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Abstract 

We review the current methods and results of lattice simulations of quantum chromodynamics at nonzero 
temperatures and densities. The review is intended to introduce the subject to interested nonspecialists and 
beginners. It includes a brief overview of lattice gauge theory, a discussion of the determination of the 
crossover temperature, the QCD phase diagram at zero and nonzero densities, the equation of state, some 
in-medium properties of hadrons including charmonium, and some plasma transport coefficients. 
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1 



I. INTRODUCTION 



Quantum chromodynamics is the well-established theory of interacting quarks and gluons. Al- 
though its Lagrangian is simple and elegant, except for high energy processes where perturbation 
theory is applicable, it is very difficult to solve. Over the past three decades ab initio numerical 
and computational methods have been devised for obtaining nonperturbative solutions. They have 
become refined to the point that a few dozen calculated quantities (decay constants, mass split- 
tings, etc.) agree with known experimental values to a precision of a couple percent [1]. These 
successes provide the opportunity to push the calculations with some confidence into new regimes 
that have not been thoroughly explored experimentally. In this review we will be interested in 
numerical simulations of strongly interacting matter under the extreme conditions of high temper- 
atures and/or high baryon number densities. 

Shortly after the big bang the universe was very likely dominated by a high temperature plasma 
of quarks, antiquarks, and gluons. As the universe expanded and cooled, hadrons emerged that 
make up today's universe. Knowing the characteristics of the plasma and the nature of the tran- 
sition to hadrons is clearly important for understanding these stages in the development of the 
universe. In the cores of some dense stars it is conceivable that the baryon number density is suf- 
ficiently high that hadrons lose their identities and merge into a plasma of quarks and gluons. The 
equation of state of such a dense plasma, for example, is important for understanding conditions 
leading to a collapse to a black hole. In heavy ion collisions at RHIC, FAIR, and soon at the LHC 
we seek to produce a quark-gluon plasma and study its properties. Since so little is known about 
the plasma, we turn to numerical simulation of high-temperature and moderate-density QCD to 
predict its properties and to guide the experimental investigation. Apart from the phenomeno- 
logical interest in such simulations, there is also intrinsic theoretical interest in understanding the 
behavior of confining field theories under extreme conditions. In particular, there are tantalizing 
predictions of still new states of matter at very high densities [2]. Lattice QCD thermodynamics is 
understandably a popular and vigorous field of research. 

Certainly, present day lattice simulations can't answer all of our questions. The current stan- 
dard methodology assumes thermal equilibrium. Moreover, simulations at nonzero densities are 
still in their infancy, so much of what we know is restricted to zero or very small baryon number 
densities. To apply lattice results to the phenomenology of heavy ion collisions requires an inter- 
mediate model, such as hydrodynamics, which takes input from lattice simulations, adds model 
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assumptions, and makes predictions about the rapidly evolving, emerging matter. For this purpose 
the most important quantities obtained from lattice simulations are the phase diagram as a func- 
tion of temperature and baryon number density, the equation of state, speed of sound, and transport 
properties, such as the viscosities and thermal conductivity. 

In this review, intended for nonspecialists and beginners, we give a brief overview of the lattice 
methodology and discuss a variety of numerical results. We discuss challenges and potential 
sources of systematic error. In Sec. |II] we give a brief introduction to lattice gauge theory and 
discuss the advantages and disadvantages of various fermion formulations. We discuss a variety 
of observables used to determine the transition temperature Tc in Sec. Inland comment on some 
disparate results. In Sec. |IV]we review our current understanding of the phase diagram at zero 
baryon density, and in Sec. |V] we do the same for nonzero baryon number densities. We discuss 
the variety of methods in current use for simulating at nonzero densities. We review the equation 
of state in Sec. |VIl In Sec. IVIII we discuss some properties of hadrons in the high temperature 
medium, and in Sec. IVIIll some results for transport coefficients. Finally, in Sec. |lX]we summarize 
briefly the current state of the field, list outstanding problems, and list some prospects for resolving 
them. 

II. THERMODYNAMICS IN LATTICE GAUGE THEORY 
A. Quantum partition function 

Quantum thermodynamics at a fixed, large volume is based on the partition function in the 
quantum canonical ensemble 



where H is the quantum Hamiltonian operator, T is the temperature, and the trace is taken over the 
physical Hilbert space. At nonzero densities the grand canonical ensemble is appropriate: 



where ^u,- is the chemical potential for the /th species and A'^, is the corresponding conserved flavor 
number. For example, in QCD we may introduce a separate chemical potential for each quark 
flavor. Zero chemical potential for a given flavor implies equal numbers of quarks and antiquarks 
of that flavor, so zero baryon number density, zero strangeness, etc. 




(1) 




(2) 
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The expectation value of an observable at temperature T is computed with respect to this 
ensemble through 



(0)=Tr 



Ocxpi-H/T + Y,IJiNi/T 



/Z. (3) 



B. Feynman path integral partition function 

The Feynman path integral formalism provides a practical basis for the computation of thermo- 
dynamic quantities, especially in quantum field theory, where there are many degrees of freedom. 
It converts the trace over quantum states into a multidimensional integration over classical vari- 
ables yfl. It is beyond the scope of this review to give a detailed derivation of the path integral 



formulation, particularly for a gauge theory with fermion fields. There are standard references 



mm. 

The classical variables in the Feynman path integral are the path "histories" of the funda- 
mental fields in Euclidean (imaginary) time x. (Imaginary, because the Boltzmann weight fac- 
tor exp(— ///r) is, in effect, a time evolution operator exp{—iHt) for an imaginary time —i/T.) 
For computational purposes the histories are discretized in x. The quantum fields at any given 
time are also discretized in three-dimensional coordinate space x. The resulting path integral is 
then a multidimensional integral over variables defined on a four-dimensional space-time lattice 
{x = X, x) . The discretization of space and time introduces an error, but the error vanishes as the 
lattice spacing is taken to zero (continuum limit). 



C. Scalar field example 

For a concrete example, consider a scalar field theory described by the Lagrange density 




where V describes the mass and self-interaction. On a hypercubic lattice with point separation a 
and a central-difference discretization of the derivatives, we can write a lattice approximation 

mx)] = ^Y^[(^{x + ali)-(^{x-ali)f + V[(^{x)]. (5) 

where /it is a unit coordinate vector in the ^ direction. A Euclidean time history is then specified 
simply by giving the values of the field (|)(x,x) on all the lattice points (x,x). Each such history 
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corresponds to a classical Euclidean action which is computed by summing its Lagrange 
density over the lattice points 

5((^)=£l[(^(x,x)]. (6) 

X,X 

The partition function then becomes a multidimensional integral over the values of the field 
(|)(x,x) at each point, weighted by the exponential of the classical Euclidean action: 

Z= fY[d<^{x,x)cxp[-Sm. (7) 

X,1 

Two important conditions on the Euclidean time history are inherited from the definition 
[Eq. ([T])] of the partition function: First, the time history x ranges over a finite interval from 
to a{Ni — 1) where 

l/T^aN^, (8) 

which establishes the relation between the temperature and the Euclidean time extent of the lattice. 
Second, to reproduce the trace over quantum states, the bosonic field {|) must be periodic under 
X ^ X + aNi. 

Similarly, the expectation value of an operator ((])), which depends on the field (|), is given by 

(0) = / ]\d^Mo{(sf)^x^[-sm/^, (9) 

where we replace the field operator ^ with its classical value when we insert O ((|)) in the integrand. 



D. QCD on the Lattice 

For a renormalizable, asymptotically free theory, such as QCD, the lattice formulation takes on 
a larger significance than just a convenient computational device. The lattice regulates ultraviolet 
divergences. The lattice constant a provides an upper bound or "cutoff" scale Tl/a for momenta. 
From this point of view the lattice formulation of the theory is every bit as respectable as other 
regularization schemes. Of course, as usual, we define the theory in the limit in which the cutoff is 
removed, i.e., a ^ 0. Before this is done all quantities we calculate have cutoff errors that vanish 
in the continuum limit. 

With the lattice regulator we apply the usual renormalization process: We select a few exper- 
imental values and use them to fix the bare parameters of the theory (quark masses and gauge 
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coupling). In this way the bare parameters depend on the cutoff (lattice spacing), but the physi- 
cal predictions should approach a cutoff-independent value in the limit of zero lattice spacing. In 
principle all regularization schemes should agree in the limit that their cutoffs are taken away. 

For QCD the fields are fermions and gluons. The groundwork for the lattice formulation of 
QCD with fermions was laid down by Wilson [7] in 1974, although there was other seminal work 
on lattice theories with local gauge invariance by Wegner Isfl, Smit, and Polyakov [|9D. To preserve 
gauge invariance, gluon variables are introduced as SU (3) matrices on the links between nearest 
neighbors of the lattice. There are four forward links per site, corresponding to the four compo- 
nents of the color vector potential A^^{x). The matrix for the link joining the site x with the site 
x + aju is then 

U^{x)=exp[igaA';,{x)Xy2], (10) 

where X^' are the Gell-Mann generators of SU (3). 

For the pure Yang-Mills theory of gluons a simple lattice form of the classical action is con- 
structed from the plaquettes Up^^y{x), i.e., the product of the link matrices surrounding the unit 
square in the forward /jv direction at site x. The single-plaquette Wilson action is simply the sum 
over all such plaquettes: 

Sg{U) = £ ^ReTr[l - C/p,„v W], (11) 
The gauge coupling = /A% appears in the coefficient 

p = 6/g2. (12) 

In the continuum limit the plaquette reduces to the familiar square of the field strength tensor 
summed over eight colors c: 

ReTr[l-[/p,^vW] - ^^Y,{F;^.f + 0{a^). (13) 

c 

In fact any closed planar loop, normalized by the area in lattice units, has the same continuum limit, 
but with a different (a^) cutoff error. For example a 2 x 1 rectangular version of the plaquette 
could also be used. If the two components are combined with the proper choice of coefficients, 
one can construct an improved gluon action that eliminates the leading cutoff correction, leaving 
errors at the next order o(a^). Relative to the leading continuum contribution, which carries the 
volume factor a^, such actions are called "tree-level (a^)" improved. Further improvements can 
even eliminate quantum cutoff corrections of the type 0{a^as). The "tadpole Liischer-Weisz" 

6 



actions [fioi [llll are in this category. Improving actions in this wa y is desirable, since it brings a 



calculation closer to the continuum limit at a given lattice spacing 111211 . 

The quark fields one for each flavor, have values on each lattice site. Since they are 

fermions, they require special treatment in the functional integration: their classical values are 
anticommuting Grassmann numbers. The fermion contribution to the action for each flavor can be 
written generally as 

SF{U,^\f) = '£W)M{U,x,yMy), (14) 

where M{U,x,y) is the Dirac matrix - essentially a lattice rendering of the familiar Dirac operator 
Ip + m. The functional integral for the partition function then becomes 

Z = J[dU][d\\f][d^exp[-SG{U)-SF{U,\\r)]. (15) 

Since the dependence on the quark fields is simply bilinear, and computing numerically with 
anticommuting numbers is nontrivial, it is standard to integrate out the quark fields immediately, 
following the rules of Grassmann integration, leaving only an integration over the gauge fields, 
weighted by the determinant of the Dirac matrix: 

Z = J[dU]exp[-SG{U)]det[M{U)]. (16) 

There are many ways to formulate a lattice fermion action, each with its advantages and disad- 
vantages. A great deal of effort over the past couple of decades has been devoted to improving the 
lattice treatment of fermions. We sketch the formulations here. For more detail, see [l4|,l5l|6|]. 

1. Wilson fermions 

The original Wilson rendering of the Dirac operator D^^,y^^, + m starts from a simple central 
difference approximation to the derivative: 

Vy\\f{x) = ■^[U^;{x)\\f{x + va) -U^{x-Va)\^{x-va)], (17) 

where the link matrices Uy{x) provide the gauge-covariance. The action constructed from this 
operator is 

^F,naive = (VyYv + m) ^(jc) . (18) 

x,V 
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It describes sixteen degenerate particles where only one is desired. Wilson remedied this undesir- 
able "fermion doubling" problem by adding an irrelevant term to the action 

5'F,Wilson = 5'F,naive - y J^Xi/(x)Av\|/(jc) , (19) 

where r is usually set to 1 and Ay,\\f{x) is the covariant Laplacian, 



Ay\\r{x) = \ [Uy{x)\\r{x + \a) + {x -\a)\\r{x -va) -2\\r{x)]. (20) 



The added term gives fifteen of the doublers masses of order of the cutoff scale I /a, leaving only 
one light state. The unwanted doublers thus become inaccessibly heavy in the continuum limit. 

It is customary to rearrange the terms in the Wilson action and multiply the field '^{x) by a 
constant to give 

^F,wiison = +Kj^\j/(x) (1 +yy)Uy{x)\\f{x + va) + (1 -yy)U^{x-va)\\f{x-va) 

x,V 

(21) 

where the "hopping parameter" K = 1/(8 + 2ma) controls the quark mass. Improvements to the 
Wilson formalism include removing tree-level O (a) errors by introducing a "clover" term in the 



action [[ibI and, for two flavors, introducing a "twisted mass" [14, isj . 

For thermodynamics applications the chief drawback of Wilson fermions has been (1) an ex- 
plicit breaking of chiral symmetry at nonzero lattice spacing, (2) a difficulty reaching low quark 
masses, and (3) a relatively poor representation of the quark dispersion relation. None of these 
difficulties is insurmountable. Chiral symmetry is restored in the continuum limit. 

It is necessary to search for the value K = K^- where the pion mass vanishes. Since this value 
depends on the inverse gauge coupling P, one gets a curve Kc((3) in bare parameter k — (3 space as 
shown in Fig. [U Lines of constant pion mass form a family of such curves (not shown) with the 
pion mass increasing as K decreases. Also shown is a high temperature crossover line Kf(|3). Its 
location depends on A^x. Where it intersects the line, we expect a true chiral phase transition. 
Pushing to stronger coupling (smaller (3) or negative quark masses (higher k) from there takes us 
into the realm of lattice artifacts: the theory has a parity broken phase at unphysical values of the 
bare parameters, as indicated. 
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FIG. 1: The bare-parameter phase diagram for two flavors of clover-improved Wilson fermions and an 
improved gauge action for zero and nonzero temperatures, illustrating the mapping necessary for thermo- 
dynamics studies with Wilson fermions. In this plot the hopping parameter K is denoted by K. The line 
of chiral critical hopping parameters (r = 0) was determined from the vanishing of the pion mass. The 
line K, indicates the high temperature crossover at A^t; = 4. It was determined from the Polyakov line (see 
section IIII Ak The region "chiral phase transition" shows where the thermal crossover happens for small 
pion masses. The parity broken phases come from lattice artifacts of Wilson fermions. The data are from 



the CP-PACS collaboration, fl, as shown in llTll . 
2. Staggered fermions 

The staggered fermion approach starts from the naive action in Eq. (fTSl) . Through a field 
transformation, the Dirac matrices can be diagonalized exactly giving four identical actions, each 
of them with one spin per site. If we keep only one of the actions, we reduce the lattice fermion 
degrees of freedom by a factor of four, which still leaves us with four fermion doublers. These 
residual degrees of freedom are called "tastes." Without further intervention, they would overcount 
the sea quark effects by a factor of four. To get approximately the correct counting, we replace the 
fermion determinant by its fourth root for each of the desired flavors: 



-stagg 



[dU] exp[-SG{U)] Yldet[Mi{U)] 



1/4 



(22) 
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In the continuum limit at nonzero quark masses, the eigenvalues of the determinant cluster in 
increasingly tighter quartets as expected from fermion doubling nloA . Then we have an SU{4) 
taste symmetry, so taking the fourth root is equivalent to using only one of them as a sea quark 
species. This procedure has generated considerable controversy. Although there is no rigorous 
proof that the method is valid, all indications so far are that the approximation is under control 
as long as we take the continuum limit before we take the quark masses to zero or fit data to a 
chiral model with taste symmetry breaking properly included [[19I1 . in which case the limits are 
completely under control. 

At nonzero lattice spacing the taste symmetry is broken, which introduces lattice artifacts. 
For example, mesons composed of a valence quark and antiquark come in nondegenerate taste 
multiplets of sixteen tastes. In the continuum limit they are degenerate. 



The asqtad 
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28D and p4fat3 12^ iBOn improvements of the stag- 



gered fermion formalism eliminate errors of 0{a^) in the quark dispersion relation and suppress 
taste splitting significantly. The asqtad suppression is somewhat better, presumably because it 
eliminates all tree level 0{a^) errors. The recently proposed HISQ action does still better [|3l|] . 
In Fig. |2] we compare the predicted and measured scaling of the splitting in the asqtad pion taste 
multiplet. 

Taste splitting can also be reduced simply by replacing the gauge link matrices in the action by 
smoothed gauge links - for example the Dublin "stout links" ll32ll . Unlike the asqtad approach, 
this method does not eliminate O (a^) errors systematically. Thus the free quark dispersion relation 
is still unimproved. 

Is taste-symmetry breaking really a problem for thermodynamics? It is believed to be most 
dramatic for the pion and less noticeable for more massive states [33]. One could argue that close 
to the crossover temperature and away from the critical point, so many excited states participate, as 
in the hadron resonance gas model, that pions do not matter much. But if we approach the critical 
point at fixed lattice spacing, taste splitting is likely to have a strong effect on the critical behavior: 
we may even get a chiral symmetry restoring transition in the wrong universality class. And 
certainly at quite low temperatures where pions dominate the statistical ensemble, taste splitting 
makes a difference. 

Taste- symmetry breaking also complicates the definition of the "physical" quark mass in a 
thermodynamics simulation. At zero temperature it is traditional to adjust the up and down quark 
masses so that the Goldstone pion (the lightest one in the taste multiplet) has the physical pion 
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FIG. 2: Plot showing that the lattice artifact taste splitting of pion masses vanishes as a^a^ in the continuum 
limit. The splitting is measured as the difference of the squared masses of the multiplet member and the 
Goldstone pion member. It is given in units of ri w 0.318 fm. The plot symbols distinguish the members 
of the multiplet. (The subscripts in the legend denote the Dirac-gamma-matrix-style classification of the 
pion tastes, ranging from singlet (s) and Yo to YoYs-) The line is drawn with unit log-log slope to test 
proportionality to a^a^. 

mass. This is legitimate, because we may restrict our attention to Green's functions whose external 
legs are the Goldstone pion. In a thermodynamics simulation, however, all members of the taste 
multiplet participate in the thermal ensemble. Thus it is more appropriate to tune the average 
multiplet mass, e.g., the rms pion mass to the physical pion mass. At a nonzero lattice spacing, 
the multiplet splitting may be so large, that goal is unreachable. In that case the physical point 
is reached only by reducing the lattice spacing together with the light quark mass. It is simply 
incorrect to claim a thermodynamics calculation is done at a physical pion mass when the rms 
mass is still much higher. 



3. Domain-wall fermions 

Neither the Wilson fermion formulation, including the clover-improved and twisted-mass ver- 
sion, nor the staggered fermion formulation are entirely satisfactory discretizations of fermions. 
Wilson fermions explicitly break chiral symmetry and its recovery requires a fine tuning. Stag- 
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gered fermions, while preserving a remnant of chiral symmetry, have a remaining doubling prob- 
lem, requiring the fourth-root trick, which is still somewhat controversial. 

A more sophisticated, somewhat indirect and more costly discretization of fermions goes un- 



der the name of "domain-wall fermions" and was developed by Kaplan [|34|] and by Furman and 



Shamir [13511 . Furman and Shamir's construction has become standard. An additional, fifth dimen- 
sion of length Ls is introduced and one considers 5-d Wilson fermions with no gauge links in the 
fifth direction, and the 4-d gauge links independent of the fifth coordinate, s, 

Sdw= £ £?(x,5)|£('Y^V^-^A^')\|/(x,5)-Mv(x,5)-P_\|/(.«,5-fl)-P+V(^^ 

(23) 

where P± = j{l iy^) are chiral projectors and we have set r = a = 1. The parameter M, often 
referred to as domain-wall height, is introduced here with a sign opposite that of the usual mass 
term for Wilson fermions [Eq. (fT9l)1. It needs to be chosen in the interval < M < 2. For free 
fermions the optimal choice is M = 1, while in the interacting case M should be somewhat larger. 
The fermion fields satisfy the boundary condition in the fifth direction, 

P-\\f{x, L,) = -m fP-\\f{x, 0) , P+\\f{x, - 1 ) = -mfP+\\f{x, L, - 1 ) , (24) 

where rtif is a bare quark mass. 

The domain-wall fermion action, Eq. (|23l) . has 5-d chiral zero modes ^ bound exponentially to 
the boundaries at 5 = and s = L^ — I, which are identified with the chiral modes of 4-d fermions 
as, 

/(x) = P+\\f{x,L, - 1) , q^{x) = P-\\f{x,0) , fix) = ^{x,L, - 1)P_ , q^{x) = W,0)P+ • 

(25) 

The left- and right-handed modes and q^ do not interact for ruf = when Ls 0° and the 
domain-wall action has a chiral symmetry. At finite the chiral symmetry is slightly broken. A 
popular measure of the chiral symmetry breaking is called "residual mass", m^^v. It is determined 
from the axial Ward identity applied at the midpoint between the two domain walls, as sketched in 
Fig- [21 This residual mass was expected to fall off exponentially in L^. But, due to lattice artifacts 
of Wilson fermions with large negative mass, there is a contribution to rtires that decreases only 



like l/Ls [|36l. l37|l . An example from a recent dynamical domain- wall fermion simulation pSl] is 
shown in Fig. HI Nevertheless, often Ls = O {10 — 20) is large enough to keep the chiral symmetry 
breaking negligibly small, especially at smaller lattice spacing (weaker coupling). 
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FIG. 3: Sketch, courtesy of Taku Izubuchi, of the domain-wall fermion setup. Left and right handed modes 
are exponentially bound to the left and right domain walls. The residual mass nii-es is determined from an 
axial Ward identity applied in the center slice. 




FIG. 4: Plot of the residual mass m,.es as a function of Lg showing its desired suppression with increasing 
Ls and increasing inverse gauge coupling p. Also shown are fits to an exponential fall-off plus a l/L^ 



contribution, from a recent 2+1 flavor dynamical domain-wall fermion simulation i38h . 



Domain-wall fermions, therefore, solve, or at least substantially alleviate, explicit chiral sym- 
metry breaking without a doubling problem. The price is a computational cost roughly a factor of 
Ls larger than that for Wilson type fermions. 

Early, Ni = 4- nonzero-temperature domain- wall fermion simulations suffered from large resid- 
ual mass, since the lattice spacing in the transition/crossover region is large, leading to much 
heavier quarks than desired [39]. More recent simulations with Nx = 8, still using Ls = 32 and 
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FIG. 5: Residual mass rrires for the recent nonzero-temperature simulations on A^t; = 8 lattices with = 32 



ll40n . At lower P, corresponding to lower temperatures, m,-ei increases rapidly, and is larger than the input 



light quai^k mass already in the transition region. 



even 96 are described in [140']. Even for the = 32 simulations with Nt = S the residual mass is 
uncomfortably large in the transition region, and getting worse at lower temperatures, correspond- 
ing to smaller [3 as shown in Fig. |5l since one would like rrires to be small compared with the input 
light quark masses. 



4. Overlap fermions 



Related to the domain-wall fermions of the previous subsection are the so-called overlap 
fermions developed by Narayanan and Neuberger [14 iL 14211 . They retain a complete chiral sym- 
metry without the doubling problem, albeit again at substantial additional computational cost. 

The overlap Dirac operator for massless fermions can be written as [|42|l . 



aDo^ = M [1 + Y5£ (ysDvv (-M) )] 



(26) 



where Dvi/(— M) is the usual Wilson Dirac operator with negative mass m = —M. As with domain- 
wall fermions < M < 2 should be used. For a hermitian matrix X, £{X) is the matrix sign 
function, that can be defined as 

e{X) = -4= ■ (27) 
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Using the fact that 8^(Z) = 1 it is easy to see that the Neuberger-Dirac operator satisfies the 
so-called Ginsparg- Wilson relation [|43l] . 



with R= 1/M. Equivalently, when the inverse of Dov is well defined, it satisfies 



(28) 



(29) 



Chiral symmetry, in the continuum, implies that the massless fermion propagator anticommutes 
with 75. As seen above, the massless overlap propagator violates this only by a local term that 
vanishes in the continuum limit. According to Ginsparg and Wilson this is the mildest violation 
of the continuum chiral symmetry on the lattice possible. Liischer [4^ has shown that any Dirac 
operator satisfying the Ginsparg-Wilson relation (l28l) has a modified chiral symmetry at finite 
lattice spacing, 

6v = /eY5 (^I-^d)v, SV = 'e¥(l-^^)Y5 ■ (30) 

or 

d\\f = iej5 1 - \\f = iejsxif , d\\f = i£YY5 , (31) 

with Y5 = Y5 (l - fjD) satisfying = 75 and, using the G-W relation, Eq. ^ = 1. 

So far only one exploratory study, on a 6^ x 4 lattice, of nonzero-temperature overlap fermions 
has been done n45\\ . The main difficulty and computational cost for overlap fermions comes from 
the numerical implementation of the matrix sign function, Eq. (l27l) . 



E. Cutoff effects 



In selecting a fermion formalism for a thermodynamics study, it is important to be aware of 
possible lattice artifacts (cutoff effects). There are two important categories of artifacts. One 
comes from an imperfect rendering of chiral symmetry. The other, from the free quark dispersion 
relation. 

It is obviously important to get the chiral symmetry right if we are simulating close to a chiral 
phase transition. Each action has its problems with chiral symmetry. For staggered fermions the 
taste splitting interferes. For Wilson fermions, the chiral symmetry is explicitly broken at nonzero 
lattice spacing. For these actions the obvious remedy is to reduce the lattice spacing. For domain- 
wall fermions, chiral symmetry is broken to the extent the fifth dimension is not infinite, and, 
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for overlap fermions, chiral symmetry is broken to the extent the matrix sign function is only 
approximated in numerical simulations. For the latter two chiral actions, this type of error can be 
reduced without also reducing the lattice spacing. 

At high temperatures where quarks are effectively deconfined, it would seem important to have 
a good quark dispersion relation, so, for example, we get an accurate value for the energy density 
and pressure. This artifact can be studied analytically for free fermions. Recently Hegde et al. ll46n 
looked at deviations from the expected free-fermion Stefan-Boltzmann relation for the pressure p 
as a function oil/N^ (equivalently a^) and chemical potential i^/T: 

f4 - L^2kP2kWTlT) i^-j (32) 

where Pzkil^/TiT) is a polynomial normalized so that P2k{0) = 1. The leading term Aq is the 
Stefan-Boltzmann term. The ratios of higher coefficients A2k/Ao measure the strength of the cutoff 
effects. These terms determine the ability of the action to approximate the continuum free fermion 
dispersion relation, and they are useful in comparing actions to the extent free quarks are relevant 
in an interacting plasma. Table IJ reproduces their results for a variety of actions. We see that the 
hypercube action ll47n has pleasingly small coefficients. The Naik (asqtad) and p4 (p4fat3) actions 
remove the second order term as designed, but the p4 action is better at sixth order. The standard 
(unimproved) staggered action (regardless of gauge-link smearing) does as poorly as does the 
standard (and clover-improved) Wilson actions. The overlap and domain-wall actions constructed 
from the standard Wilson kernel unfortunately inherit its poor behavior. Improving the kernel 
fermion action would help to reduce these cutoff effects. 



ni. DETERMINING THE TRANSITION TEMPERATURE 

We want to know the temperature of the transition from confined hadronic matter to a quark- 
gluon plasma for two obvious reasons: to interpret experimental data and to understand QCD as a 
field theory. If the transition is only a crossover, a likely possibility for QCD at the physical value 
of the quark masses as discussed below, and a true phase transition occurs only at unphysical values 
of the quark masses, then these two purposes diverge. A crossover temperature is imprecise, so its 
meaning could vary with the observable, but one can at least speak of a range of temperatures over 
which phenomenologically interesting changes take place, or one could choose one observable to 
identify a temperature. A true phase transition has a precise temperature defined by the singularity 
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A4/A0 




standard staggered 


248/147 


635/147 


3796/189 


Naik 





-1143/980 


-365/77 


p4 





-1143/980 


73/2079 


standai^d Wilson 


248/147 


635/147 


13351/8316 


hypercube 


-0.242381 


0.114366 


-0.0436614 


overlap/ 
domain wall 


248/147 


635/147 


3796/189 



TABLE I: Continuum limit scaling behavior of free massless quai^ks in various lattice formulations, based 
on an expansion [Eq. (l32l )1 of the pressure in powers of 1 /A^^ from 14611 . Shown are ratios of the expansion 
coefficients to the ideal, leading Stefan-Boltzmann coefficient. A small ratio indicates good scaling. 



of the partition function, and all observables capable of producing a signal should agree about the 
temperature. 

In this section we discuss a variety of observables commonly used to detect the transition. In 
the following sections we discuss what we have learned from them about the phase structure of 
QCD. 

Two observables are traditionally used to determine the temperature of the transition: the 
Polyakov loop and the chiral condensate. The Polyakov loop is a natural indicator of decon- 
finement. The chiral condensate is an indicator of chiral symmetry restoration. 



A. Polyakov loop and the free energy of color screening 

The Polyakov loop is an order parameter for a high-temperature, deconfining phase transition 
in QCD in the limit of infinite quark masses. At finite quark masses it is no longer an order 
parameter, but it is still used to locate the transition. It is built from the product of time-like 
gauge-link matrices. It is the expectation value of the color trace of that product: 

/ N,~l \ 

L(x,a,r) = (Tr n t^o(x,x) ). (33) 
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This quantity is gauge invariant because the combined boundary conditions for gluon and fermion 
fields require that gauge transformations be periodic under t t + Nx. Translational invariance 
insures that it is independent of x. It can be shown that the Polyakov loop measures the change in 
the free energy of the ensemble under the introduction of a static quark (excluding its mass). 

L{a,T)=exp[-FL{a,T)/T]. (34) 

In that sense the Polyakov loop is a useful phenomenological quantity as we now explain. When 
a static quark is introduced it must be screened so that the ensemble remains a color singlet. At 
low temperatures, screening is achieved by binding to it the lightest antiquark, forming a static- 
light meson. The free energy cost then consists of the self-energy of the static charge, the binding 
energy, and the self-energy of the light quark. In the quark plasma, color neutrality is achieved 
through a collective shift of the plasma charges, as in Debye screening in an ordinary electrical 
plasma. Aside from the self-energy of the static quark, which is the same at all temperatures, the 
additional free energy cost is small. So we expect Fg(a, T) to decrease abruptly in the transition 
from the confining regime to the plasma regime. 

The static-quark self energy diverges as 1/a in the limit of small lattice spacing, so it is conve- 
nient to remove it from the definitions of the free energy and the Polyakov loop: 

Ft{a, T) = Fstaticl^) +Fq{T) U^norm{T) = exp[-Fq{T) /T]. (35) 
Figure [6] illustrates the free energy from a recent lattice simulation. (Here and elsewhere, the 



temperature scale is given in MeV and in units of the Sommer parameter 1481], ''o ~ 0.467 fm. 
The latter is defined in terms of the potential V{r) between a heavy quark and antiquark. It is the 
distance where r^dV{r)/dr = 1.65.) The renormalized free energy behaves as expected. 

If we take the masses of all the quarks to infinity, we arrive at the pure SU (3) Yang-Mills en- 
semble, which has a first-order deconfining phase transition with zero L(a, T) at low temperatures 
and nonzero at high temperatures. The free energy is correspondingly infinite for T < and finite 
above. In this limit the Polyakov loop is a true order parameter for the transition. 
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FIG. 6: Renormalized screening free energy of a static quark (from the renormalized Polyakov loop) vs. tem- 
perature in MeV units (bottom scale) and ro units (top scale) for Nt = 6 and 8 from a HotQCD study com- 
paring p4fat3 and asqtad staggered fermion formulations i49l.l50l.l5lll. Measurements are taken along a line 
of constant physics with = 0. Im.,. The vertical bands here and in HotQCD figures below indicate a 
temperature range 185 - 195 MeV and serve to facilitate comparison. 

B. Chiral condensate 

1. Chiral symmetry restoration 

The second traditional observable is the chiral condensate. It is the order parameter for a 
high-temperature, chiral-symmetry restoring phase transition at zero up and down quark masses. 
At nonzero quark masses, it is no longer an order parameter, but, like the Polyakov loop, it is 
used as an indicator of the transition. It is defined for each quark flavor / as the derivative of the 
thermodynamic potential InZ with respect to the quark mass, 

(vK^M-W) = ^^ = ^(TrMri), (36) 

or the expectation value of the trace of the inverse of the fermion matrix. When the u and d quark 
masses both vanish, QCD has a i7(l) x SU (2) x SU (2) chiral symmetry, which is spontaneously 
broken at low temperatures to U{1) x SU (2), i.e., the familiar baryon number and isospin sym- 
metries. At high temperatures the full chiral symmetry is restored. The chiral condensate (vjAj/) is 
the order parameter of the broken symmetry. It is nonzero at low temperatures and zero at high 
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FIG. 7: Chiral condensate vs. temperature in MeV units (ro scale) for A^i; = 8 from [|52i1 using the asqtad 
fermion formulation. Measurements were taken along lines of constant physics with a range of light, de- 
generate up and down quark masses m,„/ specified in the legend as a fraction of the strange quark mass m^. 
An extrapolation to zero quark mass is also shown. 



temperatures. With only two flavors, the phase transition is expected to be second order, so the 
chiral condensate is continuous at the transition. When the quark masses are small, but nonzero, 
as they are in nature, the symmetry is explicitly broken and the chiral condensate does not vanish 
at high temperatures, but it is small. 

Figure |7] illustrates the behavior of the chiral condensate from a recent lattice simulation with 
two light quark flavors and one massive strange quark. Measurements were taken along "lines 
of constant physics", i.e., curves in the space of the bare parameters (gauge coupling and quark 
masses) along which the pion and kaon masses are approximately constant, whether or not they 
have their correct experimental values. The extrapolation to zero light quark mass appears to be 
consistent with the expected behavior of this chiral order parameter. 
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2. Chiral multiplets 



The restoration of chiral symmetry leads to symmetry multiplets in the hadronic spectrum. 
At low temperatures, where the symmetry is spontaneously broken, the spectrum consists of the 
familiar hadrons. At high temperatures, where the symmetry is restored, they may have analogs as 
resonant plasma excitations, at least not too far above the crossover temperature. They also control 
screening in the plasma in analogy with the Yukawa interaction. (See Sec. IVIIBi ) 

When the light quarks are massless, spontaneous symmetry breaking requires that the pion be 
massless. If the symmetry is restored at high temperatures, the pion, suitably defined as a state, 
acquires a mass. Of course, in nature, the light quarks are not massless, so the symmetry is only 
approximate, and the pion has a small mass at low temperatures. 

Another consequence of restoring the chiral symmetry with massless u and d quarks is that 
all hadronic states involving those quarks would fall into larger symmetry multiplets. Thus, for 
example, the three pions become degenerate with the /o, the three ao's become degenerate with 
the r|, and nucleons become degenerate with parity-partner nucleons. 

The classical QCD Lagrangian suggests a further U{\) chiral symmetry, which would conserve 
a flavor-singlet axial charge. This symmetry is broken at the quantum level. This quantum phe- 



nomenon is called the Adler-Bell-Jackiw axial anomaly [15311 . Whether the strength of the anomaly 
decreases in conjunction with the high temperature transition is an open question. 

If the anomaly also vanishes, the eight meson states listed above fall into a single degenerate 
supermultiplet. Again, if the light quarks are not precisely massless or the anomaly does not 
completely vanish, these statements are only approximate. 

Whether or not hadron-like resonances are observable in experiments, the multiplets appear, 
nonetheless, in calculations, most notably in simulations of hadronic screening. 

3. Singularities of the chiral condensate 

Although we require a numerical simulation to determine the chiral condensate, from general 
considerations we can predict some of its singularities at small quark mass m and small lattice 
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spacing a: 

f 

{^{a,mj)) ~ { 



Cl/2i^^T)^/m + Clm/a + analytic T < Tc, 
cim/a^ + c§m^/^ + analytic T = Tc, (37) 

c 1 m/a^ + analytic T > Tc. 



Knowing the behavior of the condensate, and in particular its singularities, is important for locating 
the phase transition. The m/a^ singularity is easily derived in perturbation theory from a one- 
quark-loop diagram. The y/m singularity at low temperatures arises in chiral perturbation theory 



at one-loop order. In this case the pion makes the loop. It is an infrared singularity caused by 
the vanishing of the pion mass at zero quark mass [54l]. Thus it appears only in the confined 
phase where the pion is massless. If we take T ^ before m ^ the square root singularity is 
replaced by the usual chiral log(m) . The term m^/^ is the expected critical behavior at the transition 
temperature. [For the expected 3-d 0(4) universality class, 5 = 0.56.] The RBC-Bielefeld group 



discusses evidence for the expected mass dependence [|54ll . 

In a calculation with three quarks with masses m„ = = and m^, it is convenient for 
comparing results of different calculations to eliminate the ultraviolet divergence by taking a linear 
combination of the light-quark and strange-quark chiral condensates 

nis 

A^,.(r) = Df,,(r)/Df,,(r = o). (38) 

The ratio ^ of the high-temperature and zero-temperature value also eliminates a common scalar- 
density renormalization factor Z5. This is the quantity plotted in Fig. [8] from a recent simulation. 
It shows the expected dramatic fall-off at the crossover. 



C. Other Observables 

Susceptibilities are often used as indicators of a phase transition. They measure fluctuations in 
the related observables. Since a transition or crossover is usually accompanied by fluctuations in 
an order parameter, the related susceptibilities tend to peak there. 



22 



1.0 



0.8 



0.6 



0.4 



0.2 



0.0 



0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 




Trn 



asqtad: N^=8 
6 

p4: N^=8 
6 



i . T[MeV] 

> * H^* ■ ^ « 



140 160 180 200 220 240 260 280 300 



FIG. 8: To give an indication of its variation with lattice spacing, we plot the chiral condensate difference 
ratio vs. temperature in MeV units (bottom scale) and ro units (top scale) for Ni = 6 and 8 from a HotQCD 



study. Results are given for both the p4fat3 and asqtad staggered fermion formulations 15 ill . Measurements 
are taken along a hne of constant physics with mud = 0. Im.,. 

1. Quark number susceptibility 

In the low temperature phase, fluctuations in quark number are suppressed by confinement for 
the same reason that the free energy of screening of a static quark is large there. At high temper- 
atures, fluctuations are common. There can also be cross-correlations. The relevant observable 
for a quark of flavor i is the expectation (Nf /V') for spatial volume V and total quark number Ni. 
This is the quark number susceptibility. It controls event-by-event fluctuations in the associated 
flavor in heavy ion collisions. For flavors / and j the generalized susceptibility (including cross 
correlations) is 

''■^-W>4g^. (39) 

We discuss the Taylor expansion of this observable in m in Sec. I VI El 

Figure [9] illustrates the behavior of the strange quark number susceptibility Xss- It shows an 
abrupt rise at the crossover. Because it has a relatively high signal to noise ratio, this quantity is 
often used to define the crossover temperature. 

We can transform the generalized quark number susceptibility Xij from the flavor basis to the 
basis in which the isospin /, hypercharge 7, and baryon number B are diagonal. The resulting 
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FIG. 9: Strange quark number susceptibility divided by the square of the temperature vs. temperature in 
MeV units (bottom scale) and ro units (top scale) for Nx = 6 and 8. Measurements are taken along a line 
of constant physics with wi,,^ = O-lm^. Results are from a HotQCD study comparing p4fat3 and asqtad 



staggered fermion formulations 15 ill . 



quantities are shown in Fig. [TOl The diagonal susceptibilities all show the expected abrupt rise 
at the crossover temperature. The offdiagonal susceptibility Xy.b shows a small nonzero value 
above the crossover. The positive correlation between hypercharge and baryon number at these 
temperatures can either be understood in terms of fluctuations in light quark degrees of freedom or 
in terms of persistent three-quark baryon states: Light up and down quarks have positive baryon 
number (1/3) and hypercharge (1/3) and their antiquarks have the opposite values. In both cases 
their fluctuations lead to a positive correlation. Strange quarks have positive baryon number (1/3) 
but negative hypercharge (-2/3). They would lead to a negative correlation, but because of their 
higher mass, they are less prevalent. So we are left with a net positive correlation. At higher 
temperatures the mass difference is irrelevant and the correlations cancel. Similar arguments can 
be made for three-quark baryonic states, where nonstrange baryons are more prevalent than strange 
baryons. 

In Fig.[lT]we show a computation of baryon (Xq) and isospin (%/) quark number susceptibilities 
from a recent computation with two flavors of clover-improved Wilson fermions on 16^ x 4 lattices 



[551], using a different normalization from that of Fig.[l0l The cutoff effects for this Wilson fermion 
simulation are seen to be significantly larger than with improved staggered fermions, as expected 
from Table ID 
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FIG. 10: Chiial susceptibility matrix in the /, Y, B basis, divided by the square of the temperature vs. tem- 
perature in units of the crossover temperature for = 6. Measurements are taken along a line of constant 
physics with niud = O.lm, from [^52^ . 
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FIG. 1 1 : Quark number susceptibilities with Wilson fermions onNi = A lattices along a line of constant 
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2. Chiral susceptibility 



The various chiral susceptibilities are based on the second derivative of the thermodynamic 
potential with respect to the quark masses 

T dHnZ 



Xi 



'J 



V dmidirij 



(40) 



For two equal mass light quarks u and d, the derivatives can be converted to expectation values of 
products of the inverse of the fermion Dirac matrix M for those species. The commonly reported 
susceptibilities are the "disconnected" chiral susceptibility 



Xdisc " 

the connected chiral susceptibility 

the isosinglet susceptibility 
and the isotriplet susceptibility 



V 



1\2' 



Xconn=^(TrM-2) 



X 



sing 



Xconn + '^Xdisc 5 



(41) 



(42) 



(43) 



Xtrip — 25^conn- 

(44) 

Figure[l2lgives an example of the peak in the disconnected chiral susceptibility at the crossover. 

Since the chiral susceptibility is the derivative of the chiral condensate with respect to quark 
mass, one can immediately derive its singularities from the expressions for the condensate in 
Eq. 



Xsing(«,m,r) ~ < 



Ci/2i'^^T)/{2y/m) +ci/a2 + analytic T < 
ci/a^ + {c^/d)m^/^^^ + analytic 
ci/a^ + analytic 



^ — ^Cl 

T > Z: 



(45) 



The RBC-Bielefeld group discusses numerical evidence for the expected mass dependence [|54|] . 
Trends in Fig.[T2]are consistent with these expectations. In the limit of zero quark mass, this quan- 
tity is infinite below the transition and finite above. In the continuum limit it has a temperature- 
independent ultraviolet divergence. Thus the BudapestAVuppertal group propose subtracting the 
zero temperature value, multiplying by the square of the bare quark mass, and dividing by the 
fourth power of the temperature [|32l] : 

m^AXdUa:'n,T)/T\ (46) 
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FIG. 12: Left panel: Disconnected light quark susceptibility vs. temperature in MeV (ro) units (bottom 
scale). Right panel: closeup of the peak region. Lines merely connect the points. Red circles and downward- 
pointing triangles, asqtad fermions. Blue squares and upward-pointing triangles, p4fat3. Squares and circles 
are along a line of constant physics with niud = O.lm^, and triangles, with = O-OSm^. All results are 
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where AXdisc i*^, T) = %(a, m, T) — %(<3, m, 0) . The cancels the scalar density renormalization 
factor. Of course, this quantity vanishes in the zero mass limit. 



D. Setting the temperature scale 

In order to quote dimensionful lattice results in physical units, it is necessary to determine 
the lattice spacing in physical units. The calibration must be based on a quantity that is reliably 
determined in zero-temperature lattice simulations. Recent favorites are the splitting of T levels, 
the mass of the Q. baryon, and the light meson decay constants, such as or /k- These scale 
determinations are not guaranteed to agree at nonzero lattice spacing and at unphysical values of 
the quark masses. Indeed, there can be substantial differences. For example, for the asqtad action 
with a nearly physical strange quark mass, a light quark mass one tenth as heavy, and a lattice 
spacing of approximately 0.12 fm, the /k scale gives a 15% lower temperature than the T splitting 
scale. For the same quark masses, at approximately 0.09 fm the discrepancy has decreased to 
8%, consistent with an approximately 0{a^) scaling error. Of course, for any quantity of interest, 
thermodynamic or not, if possible, we would like to choose a scale according to which that quantity 
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FIG. 13: Left panel: renormalized chiral susceptibi 



ity vs temperature (/k scale) from 115 TP . Right panel: 



interaction measure vs temperature (ro scale) from 15 ill . (See the definition of this quantity in Sec. IVIBl ) 



Note that the interaction measure peaks at about 20 MeV above the crossover. 

has only a small variation as the lattice spacing approaches zero. 

Recent results from Aoki et aL_\5^ give a rather different temperature Tc for the crossover 
than the HotQCD collaboration [15 IQ. Aoki et al. locate the peak in their renormalized chiral 
susceptibility at around 150 MeV (//c) for A/x = 8, 10, and 12. The HotQCD collaboration puts the 
crossover closer to 190 MeV (ro) for Ni = S and niud /m^ = 0. 1 . Here are possible reasons for the 
discrepancy: 

• Much of the difference comes from the different choice of scale. The Budapest- Wuppertal 
collaboration uses fx to set the scale, and the HotQCD collaboration uses the Sommer 
parameter ro, calibrated ultimately from T splittings 158(1 . The scale discrepancy alone could 
explain about 30 MeV of the difference. 

• Some of the discrepancy also comes from differences in lattice parameters. The Budapest- 
Wuppertal collaboration uses a smaller lattice spacing and lighter light quark mass. The 
HotQCD collaboration estimates an approximately 10 MeV (ro scale) downward shift in 
curves related to the equation of state in the continuum limit with physical quark masses. 
Some of that shift is visible in the right panel of Fig. \T3\ 

• Some may also come from differences in the fermion formulations. The Budapest- 
Wuppertal group use standard staggered fermions with stout gauge links. This approach 
reduces effects of taste splitting, but does not improve the quark dispersion relation as do 
the actions used by the HotQCD collaboration. We don't know whether such differences 
would result in a shift in a peak position, however. 
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Whatever the differences, no matter how one sets the scale, one expects all methods to give the 
same results for the same observable in the continuum limit at physical quark masses. So for now 
we are left guessing the result of taking that limit. Since most of the present difference apparently 
comes from a choice of scale, it would help our guessing to know which scale is more suitable 
for thermodynamic quantities. We have seen that the chiral susceptibility suffers from peculiar 
singularities that may make it less suitable for locating the crossover temperature. Still, the left 
panel of Fig. [13] suggests that it scales reasonably well in fK units. For the phenomenology of 
heavy ion collisions, quantities related more directly to deconfinement, such as the interaction 
measure (equation of state) and quark number susceptibility are important. As we can see from 
the right panel of Fig. [13] the interaction measure seems to show better (but still imperfect) scaling 
in the tq scale. (Preliminary HotQCD results for the chiral susceptibility are shown in the tq scale 
in Fig. [H]) 



IV, QCD PHASE DIAGRAM AT ZERO DENSITY 
A. General outline of the phase diagram 

At infinite quark mass QCD becomes a pure Yang-Mills theory, which has a well-studied. 



weak, first-order deconfining phase transition \5% . As the quark masses are decreased, the first 



order transition weakens further and devolves into a crossover, as indicated in Fig. [Ml which 
summarizes in qualitative terms the generally accepted phase structure at zero chemical potential 
in the flavors u, d, and s. 

Close to zero quark mass, chiral perturbation theory applies, and quite general arguments can 



be made about the qualitative nature of the phase transition 116011 . depending on the number of 
quark flavors with zero mass and depending on what happens to the anomaly at the transition. 
With a nonzero anomaly and only two quark flavors the transition certainly occurs at zero u and 
d quark masses, and it is in the 3-d 0{4) universality class, because of the 0{4) two-flavor chiral 
symmetry. If the strange quark is also massless, the chiral transition is first order, and, since first 
order transitions are not usually removed by small symmetry-breaking perturbations, it persists as 
the quark masses are increased. Eventually, at sufficiently large u, d, and s quark masses the system 
is too far from chiral and the first-order transition gives way to a second-order phase transition in 
the Ising or Z(2) universality class: Ising, since at nonzero quark masses, there is no remaining 
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chiral symmetry. In the m„ = vs plane a curve of such second order transitions separates 
the first-order regime from the crossover regime as sketched in the left panel of Fig. [TH 

The quantitative determination of the phase boundaries requires numerical simulation. What 
has emerged is that the second-order critical line occurs at quite small quark masses, where sim- 
ulations are particularly challenging and especially sensitive to cutoff effects [16 ll. 16211. The right 
panel of Fig. [14] shows recent results from de Forcrand and Philipsen based on a calculation using 
unimproved staggered fermions with A^^ = 4. 
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FIG. 14: Left panel: sketch of the phase diagram for QCD at zero baryon density in 2 + 1 flavor QCD as a 
function of the light quark masses showing regions where a high-temperature phase transition or crossover 
is expected. For a second-order phase transition, the universality class is shown. The physical point is 
plotted as a dot in the crossover region. Whether the expected tricritical strange quark mass m'™ is higher 
or lower than the physical strange quark mass mf^^^ is not yet firmly established. (Similar versions of this 



figure have appeared in the literature including 1631].) Right panel: result of an actual measurement of a 



portion of the 2nd order Z(2) phase boundary a.tNz = 4 from Ref. 16411 . The axes give bare quark masses in 
lattice units and the blue cross marks the physical point. 



B. Order of the phase transition for physical quark masses 

A key phenomenological question is whether there is a first order phase transition at the physi- 
cal value of the u, d, and s quark masses or there is merely a crossover. All present evidence points 
to a crossover at zero chemical potential for these species. A recent, thorough investigation has 
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been carried out by the Budapest-Wuppertal group ll65n . They examine the conventional signal 
of the peak height in the chiral susceptibility, which they renormalize using Eq. (|46|) . If there is 
no phase transition (i.e., only a crossover), the peak height should be asymptotically constant in 
the thermodynamic limit of an infinite lattice volume. If there is a first order phase transition, the 
height is infinite, but it is limited in a finite volume by finite-size effects. Asymptotically, it scales 
linearly with the lattice volume L^. If the transition is second order, the volume dependence is 
weaker, but the result is still infinite. The Budapest-Wuppertal group ran a simulation with con- 
ventional staggered fermions on stout links at = 4, 6, 8, and 10. They analyzed their data in two 
steps. First they extrapolated the inverse peak height to zero lattice spacing at fixed lattice aspect 
ratio LT, as shown in the upper panels of Fig. \T5\ Then they extrapolated the continuum values to 
infinite aspect ratio (thermodynamic limit). The result is compared in the lower panel of Fig. [T5] 
with predictions for a first order phase transition and a phase transition in the 3-d 0(4) universality 
class. The disagreement is a strong indication that there is no phase transition. 



C. Order of the phase transition for two massless flavors 

There is a related question of significant theoretical interest. When all quarks but the u and d 
are infinitely massive, we have a two-flavor theory, and, as we have observed above, as long as the 
chiral anomaly is not involved, we expect a critical point only at zero quark mass. Furthermore, 
since the two-flavor chiral symmetry SU{2) x SU{2) ~ 0(4), we expect the high-temperature 
deconfining critical point to be in the 3-d 0(4) universality class. 

This question has been investigated by several groups with somewhat contradictory results. 
Simulations with standard staggered quarks using A^i; = 4 lattices, with large lattice spa cing in the 



transition region, and hence potentially large lattice artifacts, as collected by T. Mendez [|66ll . show 
some deviations from (9(4) scaling as shown in Fig. \T6\ For 0(4) scaling, all data points should 
collapse to the curve in the figure. Two-flavor clover-improved Wilson fermion simulations jl^, 
on the other hand, indicate good 0(4) scaling as seen in Fig.flTl 

Since staggered fermions, at the large lattice spacings in the transition region on high- 
temperature lattices with small N^;, have quite large taste symmetry breaking, one might expect 
the transition to be in the U{1) y<U(l) — 0(2) universality class, rather than the 0(4) one. More 
importantly, Kogut and Sinclair [|67!] argue that finite volume effects on the fairly small (spatial) 
lattices used are quite large. Indeed they found good agreement with 0(2) scaling, when taking 
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FIG. 15: Results from 11651] . Upper panels: Inverse of the peak height in the renormalized disconnected chiral 
susceptibility vs. squared lattice spacing showing the extrapolation to zero lattice spacing. The lattice aspect 
ratio is varied from left to right. Lower panel: Inverse of the peak height in the renormalized disconnected 
chiral susceptibility v^. inverse aspect ratio cubed showing the extrapolation to the thermodynamic limit. 
Also shown are predictions for a first order phase transition and a second order transition in the 3-d 0(4) 
universality class. 



the finite volume effects into account as illustrated in Fig. 

In contradiction with the theoretical expectations and the above summarized numerical find- 
ings, D'Elia, Di Giacomo, and Pica found indications of a first-order transition using an unim- 



proved staggered fermion action and A^x = 4 116 811 . It is important to check this conclusion with 



a more refined action. One should conclude that at present the order of the high-temperature 
transition with two massless flavors is still an open question. 
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FIG. 16: A double lograrithmic plot showing strong deviations from 0(4) scaling in this parameter range 
for two flavors of staggered fermions using N-^ = 4 lattices, collected in 16611 . The function M is the chiral 
condensate (magnetization in the spin system), h the quark mass (external magnetic field) and P and 8 are 
critical exponents, t = 6/g^ — 6/gl plays the role of the reduced temperature. 
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FIG. 17: 0(4) scaling, with linear scale, for two flavors of Wilson fermions using Ni = 4- lattices, from il6h . 



D. The phase transition with a physical strange quark 



Suppose, instead, we hold the strange quark mass at its physical value and then decrease the u 
and d quark masses toward zero. According to the qualitative picture in the left panel of Fig. [I4l 
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FIG. 18: 0(2) scaling in a finite volume for two flavors of massless staggered fermions with an irrelevant 



four-fermion interaction, from 11671] . The curve comes from an 0(2) spin model simulation with "matched 
volume". 



depending on where the tricritical point lies, we could (1) encounter a critical point and enter a first 
order regime, or (2) we may have to go to zero quark mass to find a genuine phase transition. In the 
right panel we reproduce a result from de Forcrand and Philipsen suggesting the first alternative, 
but their results were obtained with an unimproved action at A^t = 4 for which we expect large 
cutoff effects. 

As we have mentioned, cutoff effects complicate the determination of the phase boundary at 
small quark mass. This is especially likely to be true for simulations based on unimproved stag- 
gered fermions (even improved staggered fermions are not entirely immune), since for them it is 
important to take the continuum limit before taking the small quark mass limit. Otherwise, one 
risks being misled by lattice artifacts. 
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V. QCD PHASE DIAGRAM AT NONZERO DENSITIES 



A. Phenomenology 

As the baryon number density is increased (i.e., all the flavor chemical potentials are increased 
from zero), according to traditional arguments, there is a chiral- symmetry restoring phase transi- 
tion along a line in the (^u, T) plane when the u and d quark masses are zero, as sketched in Fig. 



116911 . This tradition is founded on two notions. The first argues that asymptotic freedom and con- 
sequently deconfinement should reign at very high temperatures and high chemical potential. The 
second argues that spontaneous chiral symmetry breaking occurs at zero chemical potential be- 
cause, when fermions acquire a dynamical mass through symmetry breaking, the negative energy 
levels of the Dirac sea are lowered, lowering the vacuum energy. With a nonzero chemical po- 
tential the filled positive energy levels rise in energy, counteracting the advantage of a dynamical 



chiral mass, and consequently inhibiting spontaneous symmetry breaking 117011 . 

At zero u and d quark mass chiral symmetry is exact. If chiral symmetry is restored above a 
critical chemical potential and it is spontaneously broken below, analyticity requires a phase tran- 
sition. There are no such guarantees, however, when quark masses are not zero. Since we know 
from numerical simulation that at physical quark masses there is only a crossover at zero density, 
the critical line separating the chirally broken from the chirally restored phase must move away 
from the temperature axis as the quark masses are increased. It then terminates in a critical end- 
point (Te./je)- a crossover line then fills the gap from there to the temperature axis, as indicated 
by the dashed line in Fig. [191 A key phenomenological question is whether the critical endpoint is 
experimentally accessible. 



At still higher densities exotic phases have been 



proposed, including diquark condensates and 



color-flavor locked and superconducting phases |172|,lZ3|]. These phases are, thus far, completely 
beyond the reach of current lattice simulations. 



B. Lattice methods for nonzero densities 



To confirm or refute these traditional arguments requires numerical simulation. Unfortunately, 
simulations at nonzero chemical potential are very difficult, since standard lattice methodology 
requires that the Feynman path integrand be treated as a positive probability measure. In SU (3) 
gauge theory, the integrand becomes complex at a nonzero (real) chemical potential. This creates 
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FIG. 19: Sketch of the expected phase diagram for 2+1 flavor hadronic matter as a function of temperature 



and chemical potential at physical quark masses from 117 111 . The confined, chiral symmetry broken phase 
lies in the lower left, separated from the deconfined, chirally symmetric phase by a pseudocritical crossover 
line (dashed) and first order (solid) line of phase transitions. A critical point is indicated by a black hexagon. 
A nuclear matter phase transition occurs along a line extending from /u = /uq- At higher densities a color 
superconducting phase is proposed. 

a fermion "sign problem" analogous to the fermion sign problem in condensed matter physics in 
strongly-coupled electron systems away from half-filling. A solution to either problem would be 
beneficial to the other. 

To see why the problem arises, consider the naive fermion Dirac matrix M{U) = YvVy + m. 
The lattice version of the gauge-covariant derivative Vy is given by Eq. (flTl) . The terms in Vy in 
the action allow the quark to hop to next neighbor sites in the positive and negative v direction. 
Normally, hopping in all directions must have equal weight to preserve the discrete lattice symme- 
tries of axis interchange, parity, time reversal, and charge conjugation. The fermion determinant is 
then real because taking its complex conjugate corresponds to reversing the direction of hopping, 
which has the same weight. But a positive nonzero chemical potential promotes quark hopping 
in the positive (imaginary) time direction and suppresses it in the negative time direction. This is 
naturally implemented by changing the covariant time derivative as follows: 

Vo\\r{x) ^ [Uo{x)e"^\\f{x + 6a) - e^'Wl {x - Qa)\]f{x - Oa)] . (47) 
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If a quark hops along a worldline that wraps completely around the lattice in the imaginary time 
direction, it accumulates Nz factors of exp(fl/j), and the partition function receives a net enhance- 
ment exp{a/jNT:) = exp(^/r), the appropriate statistical weight for the addition of one quark to the 
grand canonical ensemble. A quark hopping backwards is interpreted as an antiquark, and its con- 
tribution is correspondingly suppressed, as it should be. With this imbalance the determinant is no 
longer guaranteed to be real. Instead it acquires a complex phase ^ °= /uV, roughly proportional 
to the lattice volume and the chemical potential. 

A complex determinant creates additional problems for staggered fermions. With 2-1-1 flavors 
of staggered fermions at nonzero densities, one requires the square root and fourth root of the 
fermion determinants. When the determinant is real, there is no phase ambiguity in the root. But 
when the determinant is complex, one has to choose the correct Riemann sheet. The ambiguities 
and an expensive remedy are discussed in [|74n . To be safe, one is limited to small /j and volumes. 

Over the years a number of methods have been proposed for treating a complex determinant. 



We give a brief account of the attempts. For recent reviews, see u5 



1. Reweighting the fermion determinant 

As a standard lattice Monte Carlo method, reweighting involves sampling the Feynman path 
integral according to one measure and then making adjustments to achieve the effect of simulating 



with a slightly different measure BTZ 



78D 



Let us see how this idea is applied to a simulation at nonzero chemical potentials /j,-, one for 
each flavor i. (To be precise, we are speaking of a quark-number chemical potential. The baryon- 
number potential is three times as large (jjibi = 3^,). The expectation value of an operator is 
given by 

(0)^= / [dU]0{U)cx^[-SG{U)]Y[Aa[Mi{U,Hi)]/Z{ti), (48) 
where = (^ui ,;U2, . . . ) and 

Z{n) = [ [dU] exp[-SG(t/)] ndet[M,-(f/,rt)]. (49) 

•' i 

Since we can't do importance sampling with the unsuitably complex determinant det[M(f/,^)] in 
the measure, we can try to do it with the real determinant det[M(i7,;U = 0)]. That is, we write 



{0)^={0R{U,n))J{R{U,^i))^, 
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(50) 



where R{U,iu) is the ratio of determinants that reweights the contributions to the integrand to 
compensate for the incorrect sampling measure: 



R(U,p)=det[M{U,ij)]/det[M{U,0)]. (51) 

Similarly, we can reweight to imitate a change in any of the parameters of the action including the 
quark masses and gauge coupling. The reweighting factor R is simply the ratio of the intended and 
actual measures. 

This procedure, often called the Glasgow method, is mathematically correct but numerically 
unstable. As the chemical potential moves away from zero, one is no longer doing importance 
sampling. In complex analysis this approach is similar to attempting to estimate a contour integral 
in the stationary phase approximation without going through the saddle point. The variance in the 
sampled values of the numerator and denominator in Eq. (|50l) grows exponentially as the lattice 
volume increases, i.e., in the thermodynamic limit. The inevitable breakdown is forestalled by 
keeping the shift in parameters small, so by working at small /j. 

A variant of this method uses the absolute value of the determinant for the sample weighting. 
The reweighting factor is then the phase [79]. This method has been applied only to small lattice 
volumes. 

Fodor and Katz propose reweighting simultaneously in the gauge coupling and fi [ISOn . They 



argue that one achieves a better overlap with this method. For example, one might expect that if 
one moves along the crossover line in the /j — T plane, the important integration domain might not 
change as rapidly as it would if one moves in some other direction. To stay on this line requires 
changing the gauge coupling along with the chemical potential. To locate the critical line, they 
follow Lee- Yang zeros of the partition function. (These zeros lie in the complex temperature or 
complex gauge-coupling plane. If there is a genuine phase transition, as the lattice volume is 
increased, they impinge on the real temperature axis and give rise to a singularity. If there is only 
a crossover, they stay harmlessly away from the real axis.) From this method they estimate the 
critical endpoint at T = 160(3.5) MeV and /jb = 3ij = 360(40) MeV at physical quark masses 



y a factor 



using conventional staggered fermions HSlll . This critical chemical potential is nearl 
of two smaller than an earlier estimate at higher quark masses and smaller volumes [|82D- Such 
sensitivity to the simulation parameters warrants further study. 
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2. Approximating the determinant with phase quenching 



With degenerate up and down quarks, simulating with the "phase-quenched" or absolute value 
of the determinant and ignoring the phase completely is equivalent to giving the up quark a pos- 
itive chemical potential and the down quark a negative chemical potential, so it is equivalent to 
simulating with an isospin chemical potential [|83n . This procedure is numerically tractable, but to 
draw conclusions regarding the phase diagram with the standard chemical potential requires some 
justification. Kogut and Sinclair present the case in ll84ll . See also liSSll . 



3. Simulating in the canonical ensemble 



Another approach is to simulate in the canonical ensemble of fixed quark (baryon) number 



m 
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88 



8911 ■ For simplicity, consider a single quark species. The canonical ensemble with 



quark number q is then obtained from the Fourier transform 

"271 



J^e""?*^ / [dU] exp[-SG(t/)]det[M(f/,^)]|^/ 



(52) 



The sign problem arises in the Fourier transform. As the quark number is increased for a given 
lattice volume and configuration, the Fourier component decreases rapidly and the sensitivity to os- 
cillations worsens, so that any discrete approximation to the Fourier transform develops a severely 
large variance. 

Meng et at. have recently proposed a new "winding number expansion" method that starts from 
the Fourier transform of the logarithm of the determinant, log(det[M(i7, /j)] ) = Trlog[M(i7, /j)] and 



proceeds via a Taylor expansion to generate the canonical partition function [|9(a 



911]. The method 



converges much better, but so far results are reported only for fairly large quark masses. 



4. Simulating with an imaginary chemical potential 



If we make the chemical potential purely imaginary, the fermion determinant becomes real, 
and a direct simulation is possible ll92n . To recover results at a physical, real chemical potential, 
we must do an analytic continuation. The success of such a continuation depends on knowing the 
analytic form of the observable as a function of chemical potential. We do if the chemical potential 
is small enough that a Taylor expansion is plausible. So in the end, the imaginary potential method 
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FIG. 20: From 18811 . Critical line as a function of quark chemical potential and temperature for four degen- 
erate flavors of unimproved staggered fermions at A'^^ = 4, bare quark mass am = 0.05, and for the spatial 
lattice volumes shown. Results from three methods are compared: the imaginary chemical potential ap- 
proach of J93I1 . the canonical ensemble approach of |88], and the multiparameter reweighting approach of 



1 8011 . A range of strong coupling values of the critical chemical potential /Jc(P = 0) is also indicated. 



provides essentially the same information as an explicit Taylor expansion about zero chemical po- 
tential. Figure [20] from de Forcrand and Kratochvila llSSn compares three methods for determining 
the critical line. Each result shown is based on the same unimproved Nf = A staggered fermion 
action. The methods agree reasonably well for /j/T < 1. Note that this is a four-flavor simulation 
with a first-order phase transition, unlike the 2-1-1 -flavor case of Fig.[ 



5. Taylor expansion method 



For small chemical potential, we may carry out a Taylor expansion of the required observables 



in terms of the flavor chemical potentials at zero chemical potential [|94l \95\\ . Since all Taylor co- 
efficients are evaluated at zero chemical potential, determining them is straightforward. However, 
the observables that give the coefficients are nontrivial. They involve products of various traces of 
the inverse fermion matrix. The traces are usually evaluated using stochastic methods. Further- 
more, as the order of the expansion grows, the number of required terms grows factorially. Thus it 



is rare to find calculations as high as eighth order [1961 . 
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6. Probability distribution function method 



The "probability-distribution-function" or "density-of-states" method is new and promis- 



ing [|98L l99n. It is related to the reweighting method. A recent variant by Ejiri combines reweighting 
with a Taylor expansion. To explain the method we start with a simple case, defining the "density 
of states" or "probability distribution function" of the plaquette P with the Wilson gauge action 
and arbitrary fermion action: 

w{P')= J[dU]d{P'-P{U))det[M{U,0)]exp[-SG{U)], (53) 

where d{P' — P) is the Dirac delta function. It is defined like the partition function, but at a fixed 
value of the plaquette. The expectation value for an observable O (P) that depends only on P is 
then 

{0)= J dP'w{P')0{P')/ j dP'w{P'). (54) 
At nonzero n we use reweighting to calculate the partition function: 

Z(a/) = j dPR{P,^)w{P), (55) 

where the plaquette-restricted reweighting function R{P,ij) is 

..p ._ ndU]d{P'-P{U))dot[M{U,^)] 
^ ^[dU]h{F-P{U))da[M{U,Q)X ^ ^ 

i.e., the ratio at nonzero and zero ^u. For the Wilson action, the gauge weight exp[— 5G(t/)] depends 
only on P, so it cancels between numerator and denominator in R{P,/j). The distribution function 
w{P) is still calculated at = according to (l53l) . 

The sign problem appears in the numeric evaluation of Ejiri offers a way to over- 

come it [|99t] . His method begins with a generalization of the distribution function, making it 
depend on three variables: the plaquette P, the magnitude of the ratio of determinants F{fj) ~ 
detM(^)/detM(0), and the phase = ImlogdetM(/j): 

w{P',\F'\,Q') = J [Jt/]5(/''-/'(t/))5(|F'|-|^l)S(e'-0)det[M(f/,O)]exp[-5G(f/)], (57) 

Note that the real, positive weight factor in the integrand comes from the /j = action. For any 
value of /J the partition function is then 

Z{/j)= J dPd\F\dQF{ij)w{P,\F\,B), (58) 
41 



where in place of the reweighting function R we now have simply F{iu) itself. 

The next step relies on the key assumption that the distribution function w{P', |F'|, 0') is Gaus- 
sian in 0. Ejiri argues that this is plausible, at least for large volume. A further assumption for 
rooted staggered fermions is that the effect on the phase of taking the fourth root is simply to 
replace by 0/4 in the Gaussian distribution. With these assumptions one can do the integra- 
tion directly, eliminating the sign problem. The result depends only on the width of the Gaussian, 
which must be determined numerically. Finally, to make the calculation of the ratio of determi- 
nants tractable, Ejiri expands log[detM(;u)] in a Taylor series in /j about ^ = 0. The same Taylor 
coefficients appear in an intermediate step in the Taylor expansion of the pressure or thermody- 
namic potential. Since one is expanding the action instead of the thermodynamic potential, the 
convergence properties are different — possibly more favorable. 

Applying this method to p4fat3 staggered fermions with the Wilson gauge action, a rather 
coarse lattice with A'x = 4, and a rather large quark mass, Ejiri locates the critical chemical potential 
at /j/T > 2.5, approximately. This is an interesting result, which awaits reconciliation with the 
questions raised by Golterman, Shamir, and Svetitsky concerning phase ambiguities of the fourth 
root of the staggered fermion determinant 1741]. 

Thus we see that all of the methods, save, perhaps, the probability distribution function method, 
are limited to quite small chemical potentials. 



7. Stochastic quantization method 



All of the above lattice methods for simulating at nonzero chemical potential evaluate the Feyn- 
man path integral using Monte-Carlo importance sampling, a technique that is inherently unstable 
when the path integrand is not positive definite. At nonzero chemical potential, the SU (3) fermion 
determinant is complex, and the wide variety of methods outlined above deal with the complex 
phase with limit ed su ccess. Instead of quantizing via the Feynman path integral meth od, A arts 
and Stamatescu lIlOOll have recently proposed using the stochastic quantization method lIlOl]. In 
the early days of lattice calculations, stochastic quantization through the Langevin equation lll02ll 



was, in fact, one of the competing numerica' 



field theory, and it met with mixed success jlOBIl . 



methods for nonperturbative calculations in quantum 



For purposes of this review, we give just a brief sketch of stochastic quantization. For a the- 
ory with a scalar field (|)(x) and action 5, we generate an ensemble of fields (|)(x, x) where x is a 



42 



fictitious Langevin time (analogous to molecular-dynamics or Markov-chain time in the standard 
importance sampling approach.) The ensemble satisfies the stochastic equation 

where ti(x, x) is a Gaussian random field (source), uncorrelated in x. As long as S has a well- 
defined minimum and we start with a solution near that minimum, without the random source the 
field relaxes to the classical solution where the action is stationary, i.e., the variational derivative 
55/6(|)(x) vanishes. The random source then induces "quantum fluctuations" about the classical 
solution. Quantum observables are estimated in the usual way as expectation values on the equi- 
librium ensemble. 

When the action S is complex, we get a complex solution and a complex stationary point, a 
region that is not reached with conventional importance sampling. The hope is that the solution is 
still attracted to the appropriate stationary point, i.e., the Langevin method is stable. Aarts and Sta- 
matescu have done some preliminary tests with simplified models that imitate the characteristics 



of QCD at nonzero chemical potential. Their results are promising 11104 



10511 



C. Curvature of the critical surface 

One question of considerable phenomenological importance can be addressed with simulations 
at small chemical potential. That is whether the Z(2) critical line sketched in Fig[14]moves closer 
to the physical quark masses as the chemical potential is increased or it moves farther away. If 
it moves closer, as shown in the left panel of Fig. [211 one may expect a true phase transition 
in a suitably baryon-rich environment, such as may occur in a moderately low-energy heavy-ion 
collision. If it moves away, as shown in the right panel, there would be no such expectation. De 
Forcrand and Philipsen set out to address this question using the imaginary chernical poten tial 



106L 



107 . 



108n . at 



method. Their results at A't; = 4 suggest that the critical line moves away [164 , 
least when all three quark flavors are close to having equal masses. 

VI. EQUATION OF STATE 

The equation of state gives the energy density, pressure, and/or entropy of the thermal QCD 
ensemble as a function of temperature at constant volume. All quantities are renormalized by 
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FIG. 21: Two possible alignments of the chiial critical surface at low chemical potential from 1641] . Left: 
the scenario permitting a first order phase transition at high densities and temperatures. Right: the scenario 
allowing only a crossover. 



subtracting their values at zero temperature. The subtraction eliminates an ultraviolet divergence, 
but the cancellation of this divergence makes the computation costly in the continuum limit, since 
one must compute the (a^^) divergent high-temperature and zero-temperature quantities inde- 
pendently and subtract them to get a finite result. 

There are two traditional methods for computing the equation of state and one recently intro- 
duced method. 



A. Derivative method 



The first method is based on the identity 



ainZ 



V dT 



(60) 



On the lattice the derivative with respect to temperature at fixed volume in the first identity trans- 
lates to a derivative with respect to 1/ (Nxat) at fixed a^, where at is the lattice spacing in the 
imaginary time direction and is the lattice spacing in the spatial direction. At fixed N-:, we 
differentiate with respect to at itself. 

For example, for the original Wilson plaquette gauge action of Eq. (fTTI) the explicit dependence 
on at and a^ goes as follows: 



Scias^at.g^) = 2/g^{as,at) 



(61) 
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where we have distinguished the timelike and spacelike plaquettes 

Pt{x) = £ReTr[l-t/p,-oW], (62) 

i 

Ps{x) = £ReTr[l-t/p,-Xx)]. (63) 

'<;■ 

In the gauge action above, we have indicated the dependence of the gauge coupling on the lat- 
tice constants a.y and at. That dependence is defined through a standard renormalization procedure 
for an anisotropic lattice: at a fixed ratio at /as and gauge coupling g, we compute an experi- 
mentally accessible, dimensionful quantity, such as the splitting of a quarkonium system. From 
the experimental value of the splitting, we can then determine the lattice constants in physical 
units. We repeat the procedure, varying g and at /as to get the full dependence of g on the lattice 
constants. 



So from Eq. (1151) with only the gauge action in this example, we have 1110911 (after setting 

at = as = a) 



aing2 



{SG/V) + {6/g^)T{Pt-Ps). (64) 



din at 

The partial derivative of the gauge coupling with respect to at is called the Karsch coefficient. 
It is known up to 1-loop order in lattice perturbation theory, but a nonperturbative calculation 
described above is necessary at experimentally accessible temperatures. As we indicated above, 
that calculation is rather involved. 



B. Standard integral method 

A second thermodynamic identity gives the pressure as the volume derivative of the thermody- 
namic potential, 

(65) 



p = T —InZ 



T 

By itself, this identity leads to an expression similar to the energy density above, but in this case 
we need the derivative of the gauge coupling with respect to the spatial lattice spacing as at fixed 
at- We have the same difficulty as before in requiring a nonperturbative calculation of an uncon- 
ventional quantity. 

But if we combine the two identities to form the interaction measure /, 

I = e-3p, (66) 
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then we get a total derivative of the gauge coupling with respect to a = as = at and the lattice 
thermodynamic identity 

V dma 

The isotropic derivative of the coupling with respect to the cutoff is just the commonly computed 
renormalization group beta function (3 = dg^/dlna. For the Wilson plaquette gauge action we get 

I=-T/V{d\ng^/d\na){SG). (68) 

So the lattice derivative is readily calculated in terms of the conventional plaquette observable and 
the beta function. With fermions present we require also the chiral condensate and the derivative 
of the quark masses with respect to the lattice spacing. These are also easily accessible in lattice 
calculations. 

We must bear in mind that the physical quantities require subtracting the zero-temperature 
values, so in the end we need the difference 

A/ = /(r)-/(0). (69) 

We will often drop the A in the following discussion and figures. 

Figure [221 shows the interaction measure difference obtained in a recent Nx = S calculation with 
equal-mass up and down quarks and a strange quark. The mass of the strange quark was held fixed 
at approximately its physical value, and the masses of the up and down quarks were set to a fixed 
fraction of the strange quark mass. Thus the temperature was varied roughly along parameter- 
space lines of constant physics, meaning light pseudoscalar mesons (at zero temperature) had 
approximately constant masses. 

To complete the determination of the equation of state, we need the energy density and pressure 
separately. The pressure is easily computed in the thermodynamic limit, in which InZ is simply 
proportional to the volume: 

InZ - -pV/T, (70) 

So the expression (|67] ) can also be written as 

^ Td{pV/T) 



V dlna 



(71) 



or, if we fix VT^ in the derivative, as 



i/T^^'-^n. (12) 

dmT 
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FIG. 22: Details of the dependence of the interaction measure on temperature in MeV units (bottom scale) 
and ro units (top scale) for three temperature ranges left to right: low, middle, and high, fo r N^ = 6 and 8 from 
a HotQCD study comparing p4fat3 and asqtad staggered fermion formulations iSlUllOll . Measurements in 
most cases are taken along a line of constant physics with = 0. Im.,. In the low temperature range the 
dashed and dash-dotted curves ai^e predictions of a hadron resonance gas model with different high mass 
cutoffs. The other curves in that range are spline fits to the data. In the high temperature range the dashed 
lines aie the leading order perturbative prediction for /Lt^ = 2nT and = nT. The brown line (the line 
passing through the points) is a fit to leading order perturbation theory plus a bag const ant, and the magenta 
line (the line passing mostly below the points) is an O {g^) EQCD prediction from jlllll . For a brief mention 
of EQCD, see Sec lylTAl 

We can then use the identity (TtTI) at fixed Nx to integrate with respect to In a (equivalently InT) to 
get the pressure: 

p{a)a'^-p{ao)a^ = - f^" M{a'){a')^ d\na' . (73) 

>/ln 

Here the lower endpoint of integration uq is a large lattice spacing, corresponding to a low temper- 
ature. If it is sufficiently low, we may take p(ao) = and the expression then yields the pressure 
at temperature T = 1/ {Nxo). 

The integration is carried out numerically, since the integrand is determined in a series of sim- 
ulations done at fixed lattice spacing. However, the spacing of the points can be set arbitrarily 
close as needed. The energy density is then obtained from z — I-\-3p and the entropy density from 
5 = £ + /?. 

This integral method was used to complete the construction of the equation of state with im- 
proved staggered quarks shown in Fig. [23l The same method has also been used in a study with 
two flavors of clover improved Wilson fermions 11121], as shown in Fig. 
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FIG. 23: Equation of state showing energy density and three times the pressure, both divided by the fourth 
power of the temperature vs. temperature for A^t; = 8. Measurements are taken along a line of constant 
physics with niud = O-lm^. Results are from a HotQCD study comparing p4fat3 and asqtad staggered 



fermion formulations 115 ih . The blue error bars on the pressure curve indicate the size of the error. The 



black bar shows a systematic eiTor from setting the lower limit of the pressure integration. 
C. Temperature integral method 

In the standard integral method above we fixed Nx and integrated E g. (1671) with respect to 
lattice spacing to get the pressure. The temperature integral method of iHM instead fixes the 
lattice spacing and "integrates" Eq. (|72|) over N-^ at fixed A^?. 

The advantage of working at a fixed lattice spacing (so fixed gauge coupling, quark masses, and 
Hamiltonian) is that the zero temperat ure s ubtraction is the same for all N-^, and we are assured 



of following lines of constant physics [I114I1 . With the standard integral method, to carry out the 
necessary subtraction, we need a separate zero temperature simulation for each high temperature 
point. Thus one may hope for a savings in computational effort. 

The disadvantage of the temperature integral method is that the integrand is known only at 
the discrete temperatures \/{Nxat) for integer A''i;. To decrease the sample interval at a given 
temperature, one must start with a smaller a,, which increases the cost substantially. Simulating 
on an anisotropic lattice helps. 

So far, the method has been tested on a pure Yang-Mills ensemble with the pleasing result 
shown in Fig.[ 
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FIG. 24: The pressure as function of T /Tpc, with Tpc the pseudocritical or crossover temperature for two- 
flavors of clove r imp roved Wilson fermions on 16^ x 4 lattices (filled symbols) and 16^ x 6 lattices (open 



symbols), from llll2ll . The simulation was done for a variety of rather heavy quark masses, indicated by the 
vector to pseudoscalar mass ratios nips/my. The lattice artifacts are larger than with improved staggered 
quarks, as expected from Table Jl 
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FIG. 25: Equation of state (interaction measure, energy density and pressure) for pure Yang-Mills theory, 
obtained u sing the T integral method at fixed lattice spacing = 0.097 fm and aspect ratio ac^/ax = 4 
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D. Step scaling method 



The standard integral method of Eq. (1731) has the disadvantage that it requires computing the 
difference between the high and zero temperature values of the interaction measure at each value 
of the gauge coupling (i.e., each high temperature point). At increasingly high temperature we 
get closer to the continuum limit and the matching zero temperature calculation becomes very 
expe nsive. Endrodi et al. propose a step scaling method that alleviates this problem to some degree 



[|ll6ll . Their idea is to compute the pressure at a given temperature as a series of differences: 

p(T)-p{0) = [piT)- p{T/2)] + [p{T/2) - p{TlA)] + .... (74) 

The increment 

p{T)=p{T)-p{T/2) (75) 

must be calculated at the same cutoff a to renormalize properly the ultraviolet divergence. In 
practice, this means matching a calculation at a given Ni=N with a calculation at A^i; = IN for the 
same bare action parameters. (The step factor 1/2 can be replaced by any factor less than 1.) The 
differences [p{T) — p{T /2)]/T^ are bounded from above, so the series 

[;.(r)-;,(0)]/r4 = ^(r)/r4 + lp-(r)/r\/2+2^p-(r)/r\/4 + ... (76) 

converges rapidly. 

Endrodi et al. suggest two ways to calculate p{T). One uses a modified form of Eq. (1731) . 

p{a,N^ = N)a^-p{a,N^ = 2N)a^ = - f'''\l{a' ,N^= N) -I{a',N^ = 2N)]{a')'^ d\na' . (77) 

Here, we have shown the Nx dependence explicitly. We assume that gq is large enough that the 
integration constants p{ao,Nx) are essentially zero. 
The second method uses the identity Eq. (ITOl ) to write 

p{T = l/{aN)) = p{a,Nx = N) -p{a,Nx = 2N) = [\nZ{Nx = 2N)-\nZ^{Nx = N)]/{N^N). (78) 

The rhs is the difference between the partition functions on two lattices of size A'^^ x 2N in which 
one lattice is intact and the other is split in half at the midpoint in imaginary time with periodic (or 
fermion-antiperiodic) boundary conditions applied to the two halves. To compute this difference, 
Endrodi et al. modify the action at the interface by introducing an interpolating parameter a such 
that a = 1 corresponds to the fully split lattice and a = 0, to the fully intact lattice. The simulation 
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FIG. 26: Circles: pressure from jllql for pure Yang-Mi 



predictions of EQCD perturbation tiieory (dotted line lllll 



Is theory at ult ra-high temperatures compared with 



lIZl, 



llSll ). The pressure is given in units of 



the Stefan-Boltzmann value and the temperature in uni ts of the temperature at the phase transition T^. The 



square is computed using the standard integral method ill9h . 



measures the derivative of InZ(a) with respect to a, which involves only fields at the interface 
The increment (l78l) is then computed from 

JlnZ(a) 



(79) 



There is still a strong cancellation involved in the integration over a, but it is a bit milder than 
the cancellation in the standard integral method. With their method they are able to reach such 
high temperatures that contact with perturbation theory is certainly expected, as sho wn in Fig.[26l 
For a lower temperature comparison of the 0{g^) EQCD prediction of Laine et al. [|lll|] with the 
interaction measure computed using standard methods, see Fig. [221 For a brief mention of EQCD, 
see Sec lyiTAl 



E. Equation of state at nonzero densities 

Heavy ion collisions involve interacting hadronic matter at relatively low baryon densities and 
high temperatures. At the other extreme, high baryon densities and low temperatures may occur 
in the cores of dense stars. In both cases we would like to know the equation of state. For the 
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low density environment of heavy ion collisions the Taylor series method is effective for lattice 
simulations. Unfortunately, thus far we have no reliable lattice method to simulate the conditions 
of dense stars. 

Consider the 2+1 flavor case of equal nonzero up- and down-quark chemical potentials /j„ = 
= /Jiij and a nonzero strange chemical potential /Js. The pressure can be expanded as follows: 

^= E<^™(n(^)"(^)"'. (80) 

n,m=0 

The coefficients c„,„ are evaluated at zero chemical potential 

1 1 1 a"+"nnz 



n\m\T^Vd{iJud/TYd{iJ,/Ty 



(81) 

^'ud.s=o 



CP symmetry requires that the coefficients vanish for odd n + mat zero chemical potential. 

For increasing n and m the coefficients c„,„ are increasingly complicated combinations of traces 
of the inverse of the lattice Dirac matrix. For a simple example, the lowest order mixed coefficient 

Such observables are technically difficult to compute because the trace is over all lattice sites 
as well as over colors. Usually such traces are evaluated by stochastic sampling methods. As the 
order n and m increase, not only are the traces more complicated, the required number of stochastic 
samples grows rapidly. In effect, the computational effort grows factorially in the expansion order. 

The quark number densities (nj,^) and (ris) can be found from first derivatives in the same 
expansion. For it is 



and for (ris). 



The leading terms in the expansion are 

^-cn(r)(^)+co2(r)(|). (85) 

The mixed coefficient cii(r) is nonzero (and negative) at low temperatures, because when we 
add a strange quark to the ensemble, it is screened by a light antiquark. This tendency persists at 
temperatures close to, but above the crossover. So for /Uud 7^ 0, the strange quark number density 
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FIG. 27: Energy density vs temperature for constant entropy per baryon number, from il2lh . 



is nonzero for = 0. In lieavy ion collisions the mean strange quark number density is zero, so 
we need to "tune" the strange quark chemical potential to obtain the experimental conditions. 

The quark number susceptibility matrix Xab for a,b E u,d,s is likewise found from second 
derivatives. For example, for the diagonal elements and the equivalent mixed light off-diagonal 
elements Xuu = Idd = Xud = Idu, we have 

The (diagonal) strange quark number susceptibility Xs = Xss is similarly obtained. The heavy-light 
mixed quark number susceptibility Xus = Xsu is 

Xus = =T 2- nmCnmiT) — j j . (87) 



n=l,m=l 



The interaction measure can also be expanded in this way |112QD- Once we have both pressure 
and interaction measure, we can determine the energy density and entropy density for any small 
chemical potential. As an example, we show the equation of state at constant entropy density per 
baryon number in Fig. [27l This is the equation of state appropriate to an adiabatic expansion or 
compression of hadronic matter, conditions that may obtain in a heavy ion collision. 
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VII, IN-MEDIUM PROPERTIES OF HADRONS 



A. Spatial string tension 

Despite its popular characterization as deconfined, high-temperature hadronic matter retains 
vestiges of confinement. Space-like Wilson loops still exhibit the area-law behavior associated 
with confinement. This is readily seen by considering dimensional reduction, in which for T ^ Tc, 
the short E uclidean time dimension (of extent l/T) is collapsed, leaving three spatial dimensions 
[Il22lll23ll . Since all dimensions are Euclidean, any one of them can be interpreted as Euclidean 
"time." We do a 90^^ rotation to turn one of the original spatial coordinates into the Euclidean time 
coordinate of a 2-1-1 dimensional field theory. 



The reduction of 4-d QCD to what is sometimes called "EQCD" [IllTll has these characteristics 



Quarks acquire a large 3-d mass <^{nT)~ + m^. This happens because the antiperiodic 
boundary condition in the small dimension requires a minimum momentum component TlT 
for that coordinate, which then contributes to the energy-momentum relation as an additional 
effective mass. 

The original fourth component of the color vector potential Aq is reinterpreted as a scalar 
Higgs-like field. The other three vector potential components become the usual vector po- 
tential of the 2-1-1 dimensional theory. We get a confining gauge-Higgs theory. 



• The 3-d and 4-d gauge couplings are related through ^3 = g4\'T. 

• The spatial Wilson loop of the original 4-d theory is now interpreted as the standard space- 
time-oriented Wilson loop of the 3-d theory. Because the theory is still confining in 3-d, we 
get a linearly rising potential with a string tension. 

In a recent calculation Cheng et al. compared the behavior of the spatial string tension of 
the full 4-d theory with predictions based on a perturbative connection between the four- and 
three-dimensional coupling and the numerically measured pr oport ionality between string tension 



and coupling in three-dimensional SU (3) Yang-Mills theory [|l24h . The comparison is shown in 



Fig. [281 The good agreement at temperatures as low as LSr^ is unexpected. 
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FIG. 28: Temperature divided by the square root of the spatial string tension vs. temperature in units of 
the crossover temperature Tq (lower scale) and in ro units (upper scale) for 2 + 1 flavors of p4fat3 quarks 
on lattices with Nx = 4, 6 and 8. The solid cur ve (w ith uncertainties indicated by the dashed lines) is the 



prediction of the dimensionally reduced theory 112411 



B. Screening masses 

The Yukawa potential can be thought of as a measure of the spatial correlation of a pion source 
and sink (the sources and sinks being static nucleons). The important insight here is that the 
screening mass m,i is the mass of a propagating particle. In the high temperature plasma we can 
consider similar correlations between interpolating operators of any type. These spatial correlators 
are controlled by confined states, as we indicated in Sec. IVII A[ Because we no longer have 
Lorentz invariance, the spatial screening masses are not expected to be equal to f reque ncies of 



real-time plasma excitations, but one can speculate that there may be a connection 112511 . In any 
case, they provide information about the structure of the plasma, they control the behavior of a 
variety of susceptibilities, and their degeneracy patterns provide information about the temperature 
dependence of symmetries. 

Euclidean thermal hadron propagators (correlators) are defined in the same way as they are at 
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zero temperature: 

Cab{x) = {OA{x)OBm , (88) 

where Oa{x) and Ob{x) are interpolating operators for the desired hadronic state. 

At zero temperature it is typical to project the correlator to zero spatial momentum, resulting in 
a time- slice correlator 

CABit)= Jd\CAB{t,x). (89) 
At large Euclidean time such a correlator has the asymptotic behavior 

Q5(0~ZAZfiexp(-M0, (90) 

where M is the mass of the hadron and Za and Zg are overlap constants. 

At nonzero temperatures one cannot explore the asymptotic limit because of the bound on 
Euclidean time < t < l/T, but one can define a spatial correlator by fixing one of the spatial 
coordinates and integrating over the other three, as in 

Cab{z) = j dtdxdyCAB{t,x,y,z). (91) 

(For fermions, it is necessary to include a Matsubara phase factor exp[/7ir?].) For large z the 
asymptotic behavior is 

CAfi(z) ~ZA(r)Z6(r)exp[-A/(r)z], (92) 

where ^{T) is the hadronic screening mass. At zero temperature ^{T = 0) = M. 

Even though the high-temperature plasma exhibits deconfining characteristics in its real time 
behavior, the spatial correlations remain confined, so the spectrum of spatial meson and baryon 
screening masses retains a gap characteristic of confinement even in the high-temperature plasma. 
However, since the screening mass for quarks approach nT at high temperatures, the valence- 
quark-antiquark meson screening masses approach InT and the valence-three-quark baryon 
screening masses approach 3nT . Furthermore, as chiral symmetry is approximately restored at 
high temperatures, they must exhibit the approximate degeneracies required by the chiral multi- 
plets. 

Armed with this background let us consider the temperature behavior of the screening mass 
/Jti{T) of the pion. At low temperatures the pion is a Goldstone boson, so the screening mass is 
small. Above the transition chiral symmetry is restored. So the screening mass rises above the 
transition temperature, approaching 2iiT. The transition temperature is marked by the change of 
slope. Figure |29] illustrates this behavior. 

56 



7 



03 4 



1 



5: 6 

o 



Free continuum 



A 







A 



A 



® 



m 







E 



00 



11 
S 



ud, N^=6 I — 

ud, N^=8 I — ^ 

us, N^=6 I — e^- 

ss, N.=6 I — ^ 



120 160 200 240 280 320 

T (MeV) 



360 400 



"T r- 



® 



ud, N^=6 I — 
ud, N^=8 I — ^ 
us, N^=6 I — e- 
ss, N.=6 I — ^ 



A 



Free continuum 



© 



m 



® ® 



© B 

[I] 



□ <I>%] 



120 160 200 240 280 320 

T (IVIeV) 



360 



400 



440 



440 



FIG. 29: Screening masses for the pseudoscalar channel (upper panel) and scalar c hann el (lower panel) 
vs. temperature in a dynamical 2+1 flavor simulation withp4fat3 staggered fermions ll2^ . Measurements 



were taken along lines of constant physics with niji ~ 220 MeV, niK = 500 MeV and Nx = 6 and 8 lll27ll . 



The isosinglet scalar /o (a) meson can be generated using the isosinglet chiral condensate 
V¥sing = W^u + W^d as the interpolating operator. It has a sizable mass at low temperature, but it 
joins the pion chiral multiplet at the transition temperature when the pion screening mass is quite 
small. Thus its mass must dip at the transition temperature and rise again, approaching InT. Thus 
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a dip in ^u/^ also marks the transition temperature. 

The chiral susceptibilities are related to hadron propagators in Euclidean space-time. For ex- 
ample, the isosinglet chiral susceptibility is 

%sing(7^) = J ^"^^(WsingWWsinglO)) = J dzQingiz^T) , (93) 

where Csing(z) is the scalar-isosinglet screening correlator generated by the isosinglet chiral con- 
densate. In addition to the /o, this correlator also contains a two-pion continuum contribution. So 
its asymptotic behavior has terms in exp(— /jy^z) as well as exp(— 2£'jtz) for Etz > Integration 
over z of these asymptotic terms yields contributions to the susceptibility that go as the inverse 
of the screening masses. At low temperatures the two-pion threshold is below the /o, so the two- 
pion continuum dominates the susceptibility. In the chiral limit this contribution is responsible for 
the 1 / -Jm singularity in the susceptibility. At high temperatures the pion screening mass rises, 
and the /o screening mass is approximately degenerate with it. Thus the two-pion continuum is 
expected to have a higher screening mass than the /o, and the susceptibility is finite in the chiral 
limit. Thus this susceptibility should be large at low temperatures and fall abruptly at the transition 
temperature. 



C. Charmonium 



To the extent the transition to a quark-gluon plasma is a crossover and not a genuine phase 
transition, one should not expect low temperature properties to change abruptly at the crossover 
temperature. Confined hadronic states may persist as plasma excitations at least for temperatures 
close to, but above the crossover temperature. One of the most studied examples is the 7/\|/, since 
it is readily observed experimentally, and, because of their large mass, charmed quarks are a good 
theore tical prob e. Numerical simulation suggests that the 7/\|/ persists to temperatures as high as 
l.STc il28Ul29|] . (See Sec. IVII below.) As the temperature increases beyond Tc, it is thought 



that screening of the heavy-qu ark potent ial eventually prevents the formation of a bound state and 
7/\|/ production is suppressed [Il29lll30n . 



1. Static quark/antiquark free energy 



There are two lattice methods for studying thermal effects in quarkoni um. The first, more 
model-dependent method, is based on a Born-Oppenheimer approximation lll30ll . One measures 
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the free energy of a static quark-antiquark pair as a function of separation r. The result is in- 
troduced into the Schrodinger equation as a temperature-dependent potential V (r, T) for a given 
heavy quark mass. As the temperature increases, screening effects weaken the potential, and even- 
tually it does not support a bound state for quarks of the given mass. This approximation should be 
good, provided the Born-Oppenheimer adiabatic approximation is good, i.e., as long as the plasma 
is able to relax to its equilibrium state on the time scale of the orbital motion of the quarks. 

Gauge invariance presents a subtlety in fashionable methods for extracting the free energy to 
be used as a Born-Oppenheimer potential. It is popular to distinguish between color-singlet and 
color-octet states of the static quark and antiquark. Since those states are supposed to be defined 
in terms of the colors of only the spatially separated quarks themselves, the separation is gauge 
dependent and probably not phenomenologically significant [|l3in . 

The potential method can be tested entirely in the context of a lattice calculation. One starts 
from the lattice static potential, derives the spectral function for the thermal quarkonium propa- 
gator (see the next subsection), and compares the result with a direct determination of the lattice 
spectral function. If the static approximation is correct, the results should agree. Recent attempts 
to follow this approach for Tc <T < l.5Tc fail to rep rodu ce any charmonium states in the spectral 
function nor any but the IS state of bottomonium lll32ll . So is the determination of the lattice 
spectral function unreliable, or is the static approximation unreliable for charmonium, or are both 
unreliable? 

Related attempts have been made to derive a heavy-quark potential suitable for use in the 



Schrodinger equation in real time (as oppose d to 



ogy is developed only in perturbation theory 1133, 



attice imaginary time), but so far the methodol- 
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2. Spectral density 



The second method is model independent, b ut mo re difficult. One measures the spectral func- 
tion of a thermal Green's function for the 7/\|/ lll36ll . The correlator is defined for some suitable 
local interpolating operator (xo,x) as 



C(xo,x,r) = (o(xo,x)o(0,0)), 



(94) 
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The spectral density p(co, q, T) is then obtained by inverting the Kubo formula for the partial 
Fourier transform C{xo, q, T) of the correlator: 

C(xo,q,T) = ^f d&pi(a,q,T)K{o),xo,T), (95) 
2% Jo 

where 

cosh(o(xo — l/2r) 

K{&,xo.T) = ■ ■ (96) 

smh((0/2r) 

Going from the Euclidean correlator C(jco, q, T) to the spectral density p(a), q, T) is a very difficult 
inverse problem. One would like to extract detailed information about the spectral density from 
quite limited information. Because of time-reflection symmetry, a simulation at A'x = 8 has only 
five, typically noisy, independent values. 

Possible remedies include (1) assuming a functional form for p and fitting its parameters {e.g., a 
delta function for the 7/\|/ or a Breit-Wigner shape), (2) decreasing the time interval at, allowing a 
larger Nx, and (3) adding further constraints on p, as in the maximum entropy method. We outline 
the last remedy in the next subsection. 



3. Maximum entropy method 



The maximum entropy method has been used to determine spectral functions in condensed 
matter physics for some time lll37n . It was first applied to lattice QCD by Asakawa, Hatsuda, 



and Nakamura [1138 



13911 . It is essentially a Bayesian method with a prior inspired by Occam's 



razor. One begins by defining an unremarkable default prior spectral density po(co, T). A typical 
choice would be the spectral density of a noninteracting quark-antiquark pair, or at least the density 
expected at asymptotically high frequency. One then requires that the spectral density p, inferred 
from the correlator data, should deviate only as much from po as the data seems to require. 

The method is applied in the context of a maximum-likelihood fit to the correlator data. We 
give a simplified description of the method. Starting from a parameterization of the spectral den- 
sity p((0, r), one predicts the correlator data and computes the usual chisquare %^[p] difference 
between prediction and data. One introduces a Shannon- Jay nes entropy for this p as follows 



Jco [p((o) - po((o) - p((o) ln[p((o)/po((o)]] . 



The "entropy" vanishes when p = Po and for small deviations from po, it is 

1 



Ja)[p(co)-po(a))]7po(ra)- 



(97) 



(98) 
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So the default prior maximizes the entropy. One then maximizes the likelihood exp(2[a, p]) or, 
equivalently, 2[a, p] itself: 

e[a,p]=a5[p]-x2[p]/2. (99) 

The positive weight a controls the balance between maximum entropy and minimum chisquare. 
In the "state-of-the-art" method, the mean of the best fits p is then obtained from the average: 

p = J d[p]daexp{Q[a,p]). (100) 

This is our answer for the spectral density. 

This method was u sed b y Asakawa and Hatsuda to study the fate of charmonium in the high- 



temperature medium 



128n. See also 



12911 and, more recently, [| 14011 . Their results for the J/\\f 
spectral density are shown in Fig. [30] and provided some of the first evidence that the J/\\f exists 
as a discernible plasma resonance for temperatures at least as high as 1.627^. before it "melts." 

When data are inadequate, results of the MEM method can be quite sensitive to the choice of 
the defau lt model. For example, one may obtain artifact excited state peaks. For some examples. 



see lll4lll . 



VIII. TRANSPORT COEFFICIENTS 



A. Shear and bulk viscosities 



Among the transport coefficients, the shear and bulk viscosities are essential to the hydrody- 
namical modeling of the expansion and cooling of the quark-gluon plasma in the aftermath of a 
heavy ion collision. They are obtained from correlators of the energy-momentum tensor at tem- 
perature T 

C^,v,pa{xo,X,T) = {T^y{xo,x)Tpa{0)) . (101) 

We need its spectral function p, which we obtain from its partial Fourier transform C^v,po('^05 T) 
and the Kubo formula 

poo 

C,jv,pa{xo,q,T) = d(£)p,jv,pa{(£>,q,T)K{(£),xo,T), (102) 

where A'(co,xo,r) is given by Eq. (l96l) . The shear (T|) and bulk (Q viscosities are obtained from 
the low-frequency behavior of the spectral function p(a), q, T) : 

,(r) = .ii„£l^^, «r) = = iimeM(5^. (103) 

co^O to 9 co^O to 
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FIG. 30: Sp ectra l density p((o) for the //\|/ for several temperatures shown in units of the crossover tem- 



perature Tc il28h . The ground state peak is visible up to 1.627^. These results are obtained in a quenched 
simulation. 



Computing the viscos ity ha s been a well known challenging problem since it was first attempted 
by Karsch and Wyld 11 14211 . The correlator is noisy, requiring high statistics. As with the J/'\\f 
correlator, this is a difficult inverse problem. A further complication is that the spectral function 
has a nasty T-independent, large (0, ultraviolet behavior p ~ (o'^, which tends to overwhelm the 
low-frequency contribution to C(x, x) for low xq. 



Possible remedies include (1) assuming a functional f orm f or p and fitting its parameters [|l42|l . 



(2) decreasing the time interval at, allowing a larger A^x 1114311 . and (3) adding further constraints 



62 



Meyer [1146 . 



147, 



on p, such as maximum en tropy 11 14411 . and working at small nonzero momentum !! 14511 . 



14811 has done a new high-statistics calculation in pure Yang-Mills theory 



and uses a parameterization of the spectral function in terms of an optimized basis set that folds 
in appropriate perturbative behavior at large co and then emphasizes deviations from this behavior. 
For the ratio of shear viscosity to entropy density, he finds r[/s = 0.134(33) at 1.657^ where 
perturbation theory gives 0.8, and for the ratio of bulk viscosity to entropy density, L,/ s < 0.15 
at 1.65rc and L,/ s < 0.015 at 3.2Tc. These results support the notion that the plasma is a nearly 
perfect fluid. 



6. Dilepton emission and related quantities 

The dilepton emission rate, the soft photon emissivity, and the electrical conductivity of the 
plasma are other important transport properties. They are obtained from the thermal correlator of 
the electric current 

Gem{xq,xJ) = (/^(xo,x)/^(0)>. (104) 
r°° d(i) 

GEM{xo,q,T)= —K{(0,xo,T)pEM{iO,q,T), (105) 

Jo 27C 

Again, this is a difficult inverse problem. The ultraviolet divergence is milder here than with the 
spectral function of the stress-energy tensor. In this case p ~ (0~. Otherwise, the same methods 
have been applied. 

The spectral density p£yvf((0, 0, T) determines the differential dilepton pair production rate 



doid^p 



, y2l7^(e^)r^,/ ^«(M (106) 

An example of the relationship between the ME M dete rmination of the spectral function and the 
resulting dilepton rate is given by Karsch et al. Il50 1 in Fig. [3T] These results show a strong 
enhancement over the free quark-antiquark pair contribution, at least up to three times Tc resulting 
from a vector meson resonance. The hard dilepton rate is obtained from the spectral function for 
co/r ^ 1, and there is rough agreement between perturbation theory and lattice simulation. 

As with the shear and bulk viscosity, the challenge is getting to low frequency to obtain the soft 



photon emissivity, and at zero frequency, the electrical conductivity: 

(107) 



co=0 
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FIG. 31: Relationship between (a) the vector meson spectral density p£A/((O)0, T) (shown here as 01/(00, T) 
and (b) the dilepton differential production rate dW / ddxP p at zero three-momentum, plotted as a fu nction 



of energy co in units of temperature for two temperatures above the crossover temperature Tc jlSOll . The 
solid lines represent a free quark-antiquark pair. The dash-dotted lines are lattice MEM results that show a 
peak corresponding to a vector meson resonance. Results are obtained in the quenched approximation. 

Extracting the spectral function itself is challenging enough. Extracting its derivative compounds 
the difficulty. Gupta et al. lIlSlll tried different Bayesian priors to constrain the spectral function. 



IX. OUTLOOK 

Numerical simulations have taught us much about the properties of high temperature strongly 
interacting matter. Here are highlights discussed in this article: 

• We have a fair understanding of the QCD phase diagram at nonzero temperature and zero 
or small baryon densities and nearly physical quark masses. 

• We have a phenomenologically useful determination of the equation of state. 

• We have a good understanding of the behavior of the quark number susceptibility. 

• We are beginning to understand the small mass limit of the chiral condensate and its related 
susceptibilities. 

• We know the plasma has persistent confining properties that are observable in screening 
masses and the spatial string tension. 

• We have some indications of the persistence of hadronic states as resonances in the plasma 
phase at temperatures close to and above Tc. 
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We are starting to determine plasma transport coefficients. 

We are starting to make contact with perturbation theory at high temperatures. 



There are many outstanding questions. Here are particularly pressing ones: 

• We need a more robust determination of transport coefficients. 

• We don't have a good way to simulate at moderately large or higher nonzero baryon number 
densities. 

• We don't know, yet, whether the critical point in the ^/T — T plane is experimentally acces- 
sible. 

• We don't know whether the tricritical point in the ms,mu4 plane lies above or below the 
physical m^. 

• We would like to understand better the behavior of the equation of state in the region where 
it overlaps with hadron resonance gas models. 

• It would be good to develop more confidence in our understanding of the continuum limit 
of phenomenologically important quantities. 

• It would be good to have high precision results from fermion formulations other than stag- 
gered for purposes of corroboration. 

• We would like to develop more confidence in our contact with perturbation theory at high 
temperatures. 

Work currently underway will help resolve some of these issues. At zero or small baryon 
number densities we expect progress with Wilson quark formulations, including clover-improved 
and twisted-mass. Simulations with domain-wall quarks will help test conclusions about chiral 
properties. Forthcoming simulations with highly improved staggered quarks (HISQ) will help 
reduce some of the lattice artifacts of the staggered fermion formulation, especially at temperatures 
leading up to Tc, where we suspect they are important. 

For simulations at nonzero baryon number densities we really need some new ideas. Perhaps 
stochastic quantization will help. For transport coefficients and the small-quark-mass region of the 
phase diagram, we may expect progress simply by applying more computing power. 
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Lattice QCD thermodynamics is a very active field. We expect continued strong progress in the 
years to come. 
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